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Symbols 

Symbols in parentheses are for computer- 
generated figures. 


A 

cross-sectional area, m 2 

A 

constant in Arrhenius law 
for reaction i 

^ix 

body force of species i in 
x coordinate direction, 
N/kg 

b{y 

body force of species i in 
y coordinate direction, 
N/kg 

Ci 

concentration of species i, 
kg-mole/m 3 

c 

speed of sound, m/s 

c p 

specific heat at constant 
pressure, J/kg-K 

Dij 

binary diffusion coeffi- 
cient, m 2 /s 

D 'p 

thermal diffusion coeffi- 
cient, m 2 /s 

E 

activation energy, cal/g- 
mole; total internal 
energy, J/kg 

F 

flux vector 

F 

1 n j 1 n 

expansion coefficients in 
Chebyshev series 

ft 

mass fraction of species i 

Gr 

Gibbs energy of reaction, 
J/kg- mole 

9i 

Gibbs energy of species i, 
J/kg-mole 

H 

source vector 


enthalpy of species i, 

J/kg 

h 0 

total enthalpy, J/kg 

h ° 

reference enthalpy at 
standard conditions, J/kg 

[I] 

identity matrix 

J 

Jacobian 


K 

equilibrium constant, 
source Jacobian 

k 

thermal conductivity, 
J/m-s-K 

kb 

reverse reaction rate 

k } 

forward reaction rate 

M 

third body 

M 

molecular weight, kg/kg- 
mole; Mach number 

N 

number of nodes 

A 

constant in Arrhenius law 
for reaction i 

Nr 

number of reactions 

N s 

number of species 

n r 

moles of species i 

P 

static pressure, Pa 

Po 

total pressure, Pa 

Q 

heat flux, J/m 2 -s 

R 

steady-state residual; gas 
constant, J/kg-K 

R° 

universal gas constant, 
J/kg-mole-K 

R°’ 

universal gas constant, 
cm 3 -atm / g-mole-K 

T 

static temperature, K 

T e 

effective temperature, K 

T n 

Chebyshev polynomial 

To 

total temperature, K 

t 

time, s 

At 

time step, s 

u 

dependent variable vector 

u 

(U) streamwise velocity, m/s 

u t 

streamwise diffusion 
velocity of species i, m/s 

v< 

diffusion velocity of 
species i, m/s 

V 

transverse velocity, m/s 

Vi 

transverse diffusion 
velocity of species z, m/s 

Wr 

species production rate of 
species z, kg/m 3 -s 


v 
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Xi 

mole fraction of species i 

X 

streamwise spatial vari- 
able, m 

Ax 

spatial step size, m 

y 

(Y) transverse spatial vari- 

able, m 

Ay 

transverse spatial step 
size, m 

7 

stoichiometric coefficient, 
ratio of specific heats 

6 

Kronecker delta function; 
central spatial difference 
operator 

V 

computational transverse 
coordinate 

X 

eigenvalue; second coeffi- 
cient of viscosity, kg/m-s 

p 

laminar viscosity, kg/m-s 

e 

computational streamwise 
coordinate 

p 

density, kg/m 3 

a 

normal stress, N/m 2 ; ef- 
fective collision diameter 


effective collision diame- 
ter, A 


T 

shear stress, N/m 2 


equivalence ratio 


diffusion collision integral 

V 

Laplacian operator 

Subscripts: 


c 

based on chemistry 

/ 

based on fluids 

i,i 

species indices 

m 

mixture 

R 

reactions, reference value 

s 

species 

sp 

evaluated spectrally 

X 

in x coordinate direction 

y 

in y coordinate direction 

Superscripts: 


n 

time level 

- 

mass-weighted value 


derivative, fluctuating 
quantity, reactants 


products 


VI 



Summary 

Research has been undertaken to achieve an im- 
proved understanding of physical phenomena present 
when a supersonic flow undergoes chemical reaction. 
A detailed understanding of supersonic reacting flows 
is necessary to successfully develop advanced propul- 
sion systems now planned for use late in this century 
and beyond. In order to explore such flows, a study 
was begun to create appropriate physical models for 
describing supersonic combustion and to develop ac- 
curate and efficient numerical techniques for solving 
the governing equations that result from these mod- 
els. From this work, two computer programs were 
written to study reacting flows. Both programs were 
constructed to consider the multicomponent diffusion 
and convection of important chemical species, the 
finite-rate reaction of these species, and the resulting 
interaction of the fluid mechanics and the chemistry. 
The first program employed a finite-difference scheme 
for integrating the governing equations, whereas the 
second used a hybrid Chebyshev pseudospectral tech- 
nique for improved accuracy. Both programs were 
used to study a spatially developing and reacting 
mixing layer, and the results were analyzed to draw 
conclusions regarding the structure of the evolving 
layer. 

1 * Introduction 

Research is currently underway, both in the 
United States and abroad, to develop advanced 
aerospace propulsion systems now planned for use 
late in this century and beyond. One such program 
is being carried out at the NASA Langley Research 
Center to develop a hydrogen-fueled supersonic com- 
bustion ramjet engine, also known as a scramjet, ca- 
pable of propelling a vehicle at hypersonic speeds in 
the atmosphere. One phase of this research has been 
directed toward gaining a detailed understanding of 
the complex flow field present in the engine over a 
range of flow conditions. Numerical modeling of var- 
ious regions of the engine flow field has been shown 
to be a valuable tool for gaining insight into the na- 
ture of these flows. This approach has been used in 
conjunction with an ongoing experimental program 
to develop an effective analysis capability (ref. 1). 

The flow field in a scramjet engine is governed by 
the Navier-Stokes equations coupled to a system of 
equations describing each of the species present ini- 
tially and produced by chemical reaction. The gov- 
erning equations were solved in prior analyses using 
either explicit or implicit finite-difference techniques, 
with the chemical reaction process modeled by an 
ideal (mixing controlled) reaction model. Using these 
approaches, analyses of various ramjet and scramjet 


configurations have been carried out, and trends that 
were established by experiments have been predicted 
(refs. 2 and 3). 

Chemical reaction is not mixing controlled 
throughout a scramjet combustor, however. Al- 
though chemical reaction may equilibrate in the rear- 
ward region of a well-designed combustor, chemistry 
in the forward portions of the combustor is certainly 
kinetically controlled. Finite-rate kinetics is, in fact, 
a critical issue in the design of flameholders in the en- 
gine, and this phenomenon must be considered along 
with the effects of molecular and turbulent fuel-air 
mixing to develop an accurate engine flow model. It 
is for this reason that attention has turned in the 
present work to a more basic and detailed analysis of 
chemically reacting flow fields. The long-term pur- 
pose of the present research is to develop detailed 
models for fuel-air mixing and reaction in an engine 
flow field and to develop accurate and efficient nu- 
merical methods for solving the equations governing 
reacting flow that result from these models. 

Because of computer resource limitations, how- 
ever, detailed modeling of the complete engine prob- 
lem cannot be considered at the present time. A 
more tractable problem that relaxes only the com- 
plexities introduced by engine geometry is posed by 
the spatially developing, primarily supersonic, chem- 
ically reacting two-dimensional mixing layer. A ma- 
jor portion of the chemical reaction taking place in a 
supersonic combustor occurs in mixing layers. All the 
difficulties introduced by the fluid mechanics, com- 
bustion chemistry, and interactions between these 
phenomena are retained by the reacting mixing layer, 
making it an ideal problem for the detailed study of 
supersonic reacting flow. 

Prior studies on supersonic reacting mixing lay- 
ers have been quite limited. A fair amount of work 
has been carried out, however, on nonreacting mixing 
layers, both supersonic and subsonic. Even without 
combustion, the results of these studies provided a 
significant amount of useful information for under- 
standing reacting layers. Carpenter (ref. 4) stud- 
ied the development of a laminar, free-shear layer 
behind steps and blunt bodies over a Mach num- 
ber range of 0 to 10. He concluded that the de- 
velopment of the layer could best be understood 
in terms of vorticity transfer. The effect of com- 
pressibility was to increase the diffusion process in 
the layer, leading to more rapid development toward 
asymptotic conditions with increasing Mach number. 
Brown and Roshko (ref. 5) studied the subsonic mix- 
ing layer that developed between nitrogen and helium 
streams and found that the layer was dominated by 
large-scale coherent vortical structures. They found 
that these structures tended to convect at a nearly 



constant speed and that the size of the structures 
and the space between them changed discontinuously 
with movement downstream by the joining of those 
structures with their neighbors. Results of their ex- 
periment “suggested that turbulent mixing and en- 
trainment was a process of entanglement on the scale 
of the large structures.” They also found that very 
large changes in the density ratio (up to 49) mea- 
sured transversely across the mixing layer had only 
a small effect on spreading of the layer. The authors 
concluded, therefore, that the significant reduction 
in supersonic mixing layer growth rate with increas- 
ing Mach number was due primarily to compressibil- 
ity effects, rather than density effects as had been 
thought in the past. 

The role of coherent structures in turbulent pro- 
cesses in mixing layers was studied further by Roshko 
(ref. 6). He found that the size of the coherent struc- 
tures and the spacing between them increased with 
increasing downstream distance. The vortices were 
found to travel at a constant speed of (i*i + U 2 )/ 2 , 
where u\ and U2 are the free-stream velocities of 
the two streams making up the layer. Each vor- 
tex also had a finite life span that began and ended 
abruptly. Coincident with two or more of these end- 
ings, a new lifespan began, with two or more vortices 
coalescing to form a new larger vortex. As noted 
above, each of these vortices was observed to move 
at a nearly constant speed, resulting in a fairly con- 
stant spacing between a vortex and its neighbor as 
they moved downstream during their lifetime. Devel- 
oped mixing layers are self-similar, however, requir- 
ing that the spacing between vortices should increase 
linearly in the mean with increasing downstream dis- 
tance. Roshko resolved this contradiction by reason- 
ing that changes in the layer must occur discontin- 
uously and irregularly along the layer such that the 
scale of the structure grew smoothly and linearly in 
the mean. Roshko further found that in the transi- 
tion region of the layer, there was only one spacing 
distance between neighboring vortices, and this spac- 
ing represented the most stable wavelength selected 
by the laminar portion of the layer. In this region 
the scales had not yet become dispersed, as they did 
further downstream in the turbulent regime. Also, 
three-dimensional effects had not come into play in 
the transition region. Finally, Roshko noted that 
mixing layer growth likely occurred not just due to 
vortex pairing, but also through an entrainment pro- 
cess by each vortex that occurred near or during 
the pairing event. Entrainment brought together 
“pieces” of fluid from either side of the layer, also 
enhancing the mixing process. Between each of these 
pairing/entrainment events, the vortices appeared to 
convect in an apparently passive fashion. 


Ferziger and McMillan (ref. 7) in studies of the 
structure in turbulent shear flows also noted the pres- 
ence of coherent structures and pairing in a devel- 
oping mixing layer. They went on to discuss the 
importance of a tearing mechanism where vortices 
tended to be torn apart by shearing and then redis- 
tributed in parts to their neighboring vortices. They 
also pointed out the importance of three-dimensional 
effects in destabilizing the layer. The coherent struc- 
tures present in the mixing layer tended to be un- 
stable to three-dimensional perturbations that de- 
stroyed the spanwise coherence of the structures. Fi- 
nally, the authors also noted that three-dimensional 
effects could also be introduced by streamwise vortic- 
ity produced by the stretching of vortical structures. 

There has been additional work in the literature 
describing important structures present in develop- 
ing mixing layers, but the authors have gone on to 
seek specific mechanisms leading to the production 
of the structures and their effect on the flow. Several 
of these authors have dealt particularly with mech- 
anisms associated with retardation of mixing in the 
supersonic development of layers. Oh (ref. 8) hy- 
pothesized that when the local mean Mach number 
exceeded 1, some fraction of the turbulence energy in 
the flow was generated by shocks that formed about 
the eddies (eddy shocks). These shocks were quite 
weak, differing little from Mach waves, but having fi- 
nite strength. Some of the eddies in the flow were 
decelerated by passing through these shocks, and 
the resulting disturbances produced pressure fluctu- 
ations. These fluctuations appeared to correlate well 
with velocity and density fluctuations in the flow. Fa- 
vorable correlations in fluctuations of pressure and 
velocity gradient gave rise to values of the pressure 
dilation term p f du^/dxj that acted as a source or 
sink of turbulent kinetic energy in the flow. This 
term vanished in incompressible flows and in low- 
speed mixing flows where there was a large density 
variation. The term took on larger values, however, 
in high Mach number free-mixing layers and acted 
as a turbulent kinetic energy sink when gradients of 
mean Mach number and density had the same sign. 
Therefore, Oh reasoned that the pressure dilation 
term could act to reduce the turbulent shear level 
in high Mach number mixing layers, thereby slowing 
the growth of the layer relative to the incompressible 
case. This effect agreed with the results cited earlier 
in this paper. Oh then carried out calculations by 
using these ideas that appeared to validate his hy- 
pothesis. Papamoschou and Roshko (ref. 9) also ob- 
served that the spreading rate of compressible mixing 
layers was significantly reduced over that of incom- 
pressible layers, and they attributed that difference 
to compressibility effects. They deduced from their 
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studies of large-scale structures in the layer that it 
was appropriate to define a natural coordinate sys- 
tem that moved with these structures. With this 
system, an alternative Mach number, termed the 
convective Mach number M c was defined as M c = 
(u — u c )/a, where u is the free-stream velocity, u c is 
the convection velocity of the large-scale structures, 
and a is the local speed of sound. The reduction 
in mixing layer spreading rate (by approximately a 
factor of 3 or 4) was shown to correlate well with 
increasing convective Mach number beginning with 
M c ss 0.5 and leveling off for M c > 1.0. Reduced 
spreading therefore seemed to the authors to be due 
to a stabilizing effect of the convective Mach number. 

Hussaini, Collier, and Bushnell (ref. 10) offered 
a possible explanation for the correlation of mixing- 
layer spreading rate with convective Mach number. 
Their explanation was tied to the formation of the 
eddy shocklets that were described earlier. The au- 
thors studied numerically the behavior of an eddy 
convecting subsonically, relative to a locally super- 
sonic flow, with a convective Mach number greater 
than one. Such flows could therefore support tran- 
sient shock structures associated with the eddies. 
As the eddy accelerated in the supersonic flow, an 
eddy shocklet formed which tended to distort the 
eddy. As this process continued, an eddy bifurca- 
tion occurred, resulting in the formation of a vor- 
tex of opposite circulation. Additionally, the length 
scale of the original vortex was reduced. Therefore, 
it was seen that eddy shocklets could reduce turbu- 
lent mixing through both the production of counter 
fluctuating vorticity and the reduction of turbulence 
scale. The authors stated that the mechanism for 
these effects resulted from the instantaneous inviscid 
pressure field induced about the front of the eddy. 
The authors further noted that the induced pressure 
field would always counter the initial vortex circu- 
lation over a portion of its contour, and for long 
enough times and weak enough eddies, the forma- 
tion of counter vorticity and consequent eddy split- 
ting would occur, resulting in a significant alteration 
of the mixing-layer structure. 

Many, if not all, of the important features de- 
scribed above for nonreacting subsonic and super- 
sonic mixing layers also occurred in reacting lay- 
ers. A majority of the studies on reacting mix- 
ing layers were carried out at subsonic rather than 
supersonic speeds, however. Yule, Chigier, and 
Thompson (ref. 11) found that, consistent with non- 
reacting flow, many combusting flows contained co- 
herent burning structures that interacted as they 
were convected downstream. They termed the burn- 
ing region associated and moving with an eddy a 
“flamelet” and found that the flamelet formed in 


only part of an eddy. They found a range of eddy 
types existed in a diffusion flame (that occurred in a 
nonpremixed reacting mixing layer). Initially there 
existed unstable laminar flow that contained an un- 
stable laminar diffusion flame. That region was fol- 
lowed by one containing sheets of vortex rings with 
smooth tongues of flame at the interfaces between 
the vortices and unburned reactants. This region 
was followed by a zone of other orderly vortex struc- 
tures, including helical vortices, which also produced 
relatively smooth tongues of flame. This zone con- 
tained the characteristics of transition observed in 
nonreacting flow. Here, viscous forces have a stabi- 
lizing influence on the flow. As the viscous forces 
became less important and inertial forces predomi- 
nated further downstream, the authors found that 
the orderliness of the eddies decreased and the flow 
became increasingly unstable and three-dimensional. 
With the introduction of three-dimensional effects, 
randomly moving cell-like flamelets also appeared. 
Even further downstream, this process evolved into a 
fully turbulent flow with eddies containing coherent 
ragged regions of burning, forming islands that were 
completely separated from the main flame. Yule et 
al. (ref. 11) also examined the structure of a single 
eddy containing a flamelet in a simple gas diffusion 
flame. The basic structure of a transitional eddy be- 
fore it interacted with other eddies is given in figure 1, 
which was taken from reference 11. The eddy con- 
tained separate regions of fuel and air that rolled up 
into the vortex, as well as a viscous core containing 
a mixture of fuel, oxidants, and products. A flame 
existed along the interface region where large trans- 
verse gradients of temperature and species concen- 
tration occurred. The local thickness of this region 
depended on the residence time and strength of the 
vortex, the local diffusion coefficients, and chemical 
kinetics. The molecular mixing required before fuel 
and air react was enhanced in the eddy by stretch- 
ing of the fuel/air interface due to the vorticity that 
the eddy contained. Preheating of fuel and air then 
took place primarily along the interface zone where 
mixing was taking place on a molecular scale. Com- 
bustion then occurred in the interface at or near sto- 
ichiometric conditions. During these processes, the 
vortex continued to convect downstream, and the in- 
duced velocity within the eddy due to its vorticity 
continued to produce valleys and an increase in vor- 
tex dimensions. This eddy growth resulted in further 
entrainment of fuel and air, producing flame and mix- 
ing layer growth. 

Yule et al. (ref. 11) then went on to discuss the 
evolution of turbulent eddies from transitional ed- 
dies. The structure is pictured in figure 2, again 
taken from reference 11. The eddy has now taken 


3 



on a three-dimensional structure, and it has begun 
to lose the circumferential coherence about its asso- 
ciated flamelet. Additionally, there now existed an ir- 
regular vorticity distribution within the eddy, which 
was interpreted to be due to the presence of smaller 
eddy scales now existing within the main eddy. Mix- 
ing down to molecular levels was still produced by 
vortical stretching, and the process appeared, in fact, 
to be more pronounced in the turbulent eddy. In ad- 
dition, the irregularity of the structure also produced 
a range of flamelet structures, resulting in a “ragged” 
flame front trailing the eddy. The authors concluded 
their study of large coherent structures in reacting 
flow by noting that such structures could lead to 
overall reduced combustion efficiency because of un- 
mixedness. Unmixedness occurred when fuel and air 
could not effectively mix because each gas was bound 
up in vortical structures during its passage through 
a combustion region. They did suggest that large 
eddies could be broken up by increasing the shear 
stresses in the flow in regions of steep velocity gradi- 
ent or by the imposition of swirl into the flow. 

Masutani and Bowman (ref. 12) also studied the 
structure of a chemically reacting plane mixing layer. 
They examined the reaction in the mixing layer be- 
tween a stream of dilute nitrous oxide and a stream 
of dilute ozone and observed similar behavior to that 
seen by the previous authors. They found that the 
mixing layer had three streamwise states. First, 
there existed fingers of unmixed free-stream fluid 
that sometimes reached entirely across the layer. 
Next, there appeared a region of mixed fluid in a 
finite-thickness interfacial diffusion zone that bor- 
dered parcels of unmixed fluid. Finally, the layer 
consisted of regions of mixed fluid of nearly homoge- 
neous composition in a global sense. 

Keller and Daily (ref. 13) conducted an experi- 
mental study of a gaseous, two-stream, reacting mix- 
ing layer flow fueled by propane, with one stream 
made of hot combustion products and the other 
stream containing cold unburnt reactants. They 
found that the mixing layer structure was qualita- 
tively unaffected by heat release for the range of con- 
ditions that they studied. Mungal, Dimotakis, and 
Hermanson (ref. 14) experimentally studied the re- 
acting mixing layer created between a dilute hydro- 
gen stream and a dilute fluorine stream over a wide 
range of conditions. They also observed the pres- 
ence of large hot coherent structures in the layer that 
strongly influenced the mixing and entrainment of 
fuel and oxidant and the overall structure of the flow 
field. 

Hermanson, Mungal, and Dimotakis (ref. 15) ex- 
tended the work described in reference 14, but with 
significantly higher heat release. They found that at 


the higher temperature resulting in this case, the flow 
still appeared to be dominated by large-scale struc- 
tures that were separated by cold tongues of fluid 
that extended well into the layer. Thus, the struc- 
ture did not appear to be altered by heat release 
and continued to be predominantly two-dimensional. 
They also found that with significant heating and 
the resulting large density changes, the shear layer 
thickness did not increase, and in fact showed a 
slight decrease. This reduction in layer thickness 
with increasing heat release was further confirmed 
by the resulting velocity profiles that showed no- 
ticeably higher values of transverse velocity gradi- 
ent with increased heating. The authors then went 
further to note that since the layer width did not in- 
crease with temperature, and since the density of the 
layer was substantially reduced by heating, the volu- 
metric entrainment rate of free-stream fluid into the 
layer must also be greatly reduced by heat release. 
Pitz and Daily (ref. 16) carried out an experiment 
to study a turbulent propane-air mixing layer down- 
stream of a rearward facing step. They also found 
that large-scale structures dominated the flow and 
that the growth of these eddies influenced the reac- 
tion zone. Reaction took place mainly in the eddies, 
although the eddies were not confined to the velocity 
gradient region of the layer. Therefore, the result- 
ing flame spread faster into the premixed reactants 
than did the mixing layer defined by the mean ve- 
locity. Thus, the region of the mixing layer defined 
by the velocity gradient did not coincide with the 
region of high chemical reaction and heat transfer. 
Broad well and Dimotakis (ref. 17) surveyed a number 
of recent papers describing experiments on reacting 
mixing layers. Based on these papers and their expe- 
rience, they then discussed the implications for mod- 
eling such flows. Their three principal conclusions 
were that molecular transport retained a significant 
role in turbulent mixing phenomena, even when the 
flow was fully developed; large-scale structures con- 
trolled entrainment, which then provided conditions 
for the subsequent mixing processes; and mixing lay- 
ers remained unsteady at the largest temporal and 
spatial scales. 

Reacting mixing layer studies using analytical or 
numerical approaches have also been carried out. 
Carrier, Fendell, and Marble (ref. 18) used a sin- 
gular perturbation technique to modify their Burke- 
Schumann thin flame solution for a more realis- 
tic finite-thickness reaction zone in a mixing layer. 
They studied the effect of fluid strain on the flame; 
their strain increased the interfacial exposure of fuel 
and oxidant, and convected additional reactant into 
the flame. Riley and Metcalfe (ref. 19) directly 
simulated a subsonic, temporally developing and 
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reacting mixing layer by using a pseudospectral nu- 
merical method and a binary single-step irreversible 
reaction with no heat release. Using this approach, 
they were able to consider the effect of the turbulence 
field on chemical reaction. Their results were shown 
to be consistent with similarity theory and in approx- 
imate agreement with experimental data. McMurtry, 
Jou, Riley, and Metcalfe (ref. 20) extended the pre- 
ceding work to consider the effect of chemical heat 
release on a subsonic, temporally developing mixing 
layer. They solved both the compressible form of 
the governing equations as well as a more computa- 
tionally efficient form of the equations valid for low 
Mach numbers. Reaction was again modeled with a 
binary, single-step, irreversible reaction. The authors 
found with their simulations that the thickness of the 
mixing layer and the amount of mass entrained into 
the layer decreased when the heat release rate due to 
exothermic reaction was increased. Likewise, the re- 
sulting product formation also decreased as the heat 
release rate increased. 

Menon, Anderson, and Pai (ref. 21) studied the 
stability of a laminar, premixed, spatially develop- 
ing, supersonic mixing layer undergoing chemical re- 
action. They introduced an infinitesimal disturbance 
into the layer and examined its spatial stability for 
both reacting and nonreacting flows. Chemical reac- 
tion was shown to have a significant effect on flow 
stability. The authors found that with reaction, 
the disturbance amplification rate was higher and 
the wave speed lower as compared with nonreactive 
cases. Also, the free-stream Mach number was shown 
to have little effect on stability when the flow was 
reacting. 

In this study, a numerical model has been de- 
veloped for describing general two-dimensional, high 
subsonic or supersonic, chemically reacting flows. 
This model was then adapted to a supersonic, chem- 
ically reacting mixing layer. Reaction in many prac- 
tical devices takes place in mixing layers, so that 
the problem chosen, while being geometrically sim- 
ple, still retained the fluid mechanical and chemical 
complexities that were under consideration. Com- 
puter programs have been developed that numeri- 
cally solve the governing equations resulting from 
the model. The programs used either a modi- 
fied MacCormack technique or a hybrid Chebyshev 
pseudospectral technique to solve the Navier-Stokes 
and species continuity equations that describe mul- 
tiple species undergoing chemical reaction. Momen- 
tum, heat, and mass diffusion were described by laws 
based on kinetic theory; chemistry was defined with 
a multicomponent finite-rate scheme; and a real gas 
thermodynamics model was employed. 


Using the computer programs developed in this 
work, detailed studies of the supersonic, spatially 
developing and reacting mixing layer were under- 
taken. The accuracy of the finite-difference and spec- 
tral programs was compared for a simple test case. 
No attempt was made at this point, however, to 
choose the preferred approach. Several phenomena 
observed only for subsonic reacting mixing layer flows 
were then sought in the supersonic layer by using 
both methods. The studies were undertaken, first, 
to verify the existence of the phenomena, and sec- 
ond, to explore the effect of the phenomena upon the 
development of a supersonic layer relative to that 
observed in the subsonic layer. 

Because of their importance in subsonic layers, 
consideration was given first to the existence of vor- 
tical structures in a supersonic reacting mixing layer. 
The effects of such structures on the development 
of the layer were then explored and compared with 
the literature cited earlier. Particular emphasis was 
given in this study to the mixing of fuel and oxidant 
in the layer, the resulting chemical reaction, the ef- 
fect of chemical heat release on mixing, and the ex- 
istence of supersonic unmixedness. The stability of 
a supersonic reacting mixing layer was also explored 
in this work. The existence of a transition zone in 
a particular mixing layer configuration was first con- 
sidered, and then mechanisms necessary to produce 
transition of the layer were examined. Emphasis was 
also given to the effects of transition on fuel-oxidant 
mixing and chemical reaction in the zone, and the ef- 
fects of chemical heat release were again considered 
in this region. 

The development of the numerical methods em- 
ployed in this study is given in section 2 along with 
appropriate calculations to check the methods. Sec- 
tion 3 describes extensions of the methods developed 
in section 2 to two dimensions. The detailed physical 
models used to describe the complex reacting flows 
to be studied are also described in this section. Fi- 
nally, section 4 describes studies of several supersonic 
reacting mixing layer cases using the finite-difference 
and spectral computer codes that were developed. 
Conclusions resulting from the mixing layer studies 
of section 4 are given in section 5. Based on these 
conclusions, directions for further research in super- 
sonic chemically reacting flows are then discussed. 

2. Development of Numerical Methods 

The numerical methods to be employed for study- 
ing chemically reacting flows are developed in this 
section. Two classes of algorithms are developed, 
the first based on established finite-difference tech- 
niques and the second based on spectral techniques. 
Spectral schemes are high-order methods and offer 
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the high level of accuracy required for combustion 
studies. These methods have been used quite suc- 
cessfully in studies of transition of flows from laminar 
to turbulent states, problems not unlike those to be 
considered in this work. 

To solve the equations governing chemically re- 
acting flows, the spatial derivatives must first be 
discretized, and then an appropriate temporal dis- 
cretization must be chosen in order to advance the 
equations ahead in time. The temporal scheme must 
be chosen carefully because the system of partial 
differential equations describing chemically reacting 
flows can be stiff because of the highly disparate time 
scales that exist among the equations. Certain chem- 
ical reactions in a combustion kinetics system can 
take place on an extremely short scale of the order 
of ltr 9 seconds, whereas the fluid dynamics may re- 
quire from 10“ 3 to 10 seconds for a typical case to 
reach steady-state conditions. (Time scales as small 
as 10 -12 seconds were observed in these studies, but 
these scales were later found to arise from nonphys- 
ical behavior of certain global chemistry models at 
early integration times.) There are, of course, sev- 
eral intermediate scales lying between these two ex- 
tremes. Mathematically, stiffness is often defined by 
examining the eigenvalues of the Jacobian of the gov- 
erning equation system and noting that the ratio of 
the real part of the largest to real part of the small- 
est eigenvalue is a large number. The former physical 
definition is perhaps the more useful test of stiffness; 
it is felt directly in the numerical integration of stiff 
systems through the required proper choice of the in- 
tegration time step. This requirement will be dealt 
with now, and then a discussion will follow concern- 
ing integration of the spatial part of the problem. 

Stiffness in the system of equations governing 
chemically reacting flows typically arises from the 
source terms in the equations describing production 
and loss of the chemical species that are present. 
Large values for these source terms produce rapid 
changes in the dependent variables being sought and 
result in the very short time scales discussed in the 
previous paragraph. To explore the problem of mixed 
(short and long) time scales, consider the ordinary 
differential equation (ODE) system (ref. 22) 

f = w (i » 

where f = [/i,/ 2] T , f(0) = [2, l] r , and 
'-500.5 499.5' 

A = 

499.5 -500.5 


The eigenvalues of [A] are Ai = —1000.0 and A 2 = 
— 1.0, and the solution to equation (1) follows as 

h(t) = 1.5e -t + 0.5e -looot ) 

\ ( 2 ) 

f 2 {t) = 1.5e -t - 0.5e -looot j 

Note that the solutions f\ and fa have a rapidly de- 
caying component corresponding to Aj and a much 
more slowly decaying component corresponding to 
A 2 . If this problem were solved numerically, accu- 
racy would require that the solution be advanced 
from the initial conditions by using very small time 
steps. However, once the solution dominated by Ai 
decays, it is preferable to advance the solution by 
using larger time steps that would still maintain an 
acceptable level of accuracy. Care must be taken 
in picking a numerical algorithm that will allow this 
choice of time step. Otherwise, the numerical stabil- 
ity of the solution continues to be dictated by Ai even 
though its component has decayed, and very small 
time steps are still required to maintain stability. In 
response to this difficulty, several authors, includ- 
ing Bussing and Murman (ref. 23); Stalnaker et al. 
(ref. 24); Widhopf and Victoria (ref. 25); and Smoot, 
Hecker, and Williams (ref. 26) recognized that the 
stiff source terms in the system of equations govern- 
ing chemically reacting flow should be evaluated im- 
plicitly. Therefore, for these studies, algorithms are 
developed with the source terms written implicitly 
at the new time level in the integration step. Other 
terms in the governing equations that do not lead to 
stiffness can still be evaluated explicitly (refs. 23-26). 

Next, the computation of spatial derivatives in the 
governing equations is considered. The importance 
of accurately modeling spatial derivatives cannot be 
overemphasized. Chemical reaction does not take 
place until fuel and oxidant are brought together and 
macroscopically mixed by convective transport and 
then mixed down to the microscopic (molecular) level 
by diffusive processes. To model these processes, spa- 
tial derivatives must be accurately computed. Be- 
cause of computer storage limitations, higher order 
numerical methods were indicated. 

Higher order finite- difference schemes offered one 
option for computing the spatial derivatives. An- 
other option was apparent from earlier studies where 
methods were developed for computing highly accu- 
rate solutions of the Euler equations. In one study, 
Hussaini et al. (refs. 27 and 28) used a spectral 
collocation method to compute the required spatial 
derivatives in the governing equations. With this ap- 
proach, several problems governed by the Euler equa- 
tions were successfully solved and accurate solutions 
were obtained on relatively coarse grids as compared 
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with finite-difference solutions of the same problems. 
Spectral methods are based on the representation of 
the solution to a problem / by a finite series of global 
functions X of the form 

TV 

f[x ) = ) ^ a n X n (x) (3) 

n=0 

where a n are the expansion coefficients of the se- 
ries (ref. 29) and X n should be a complete orthog- 
onal set. Spatial derivatives of / are then approxi- 
mated by taking derivatives of the corresponding se- 
ries (eq. (3)). If properly applied, the high order 
approximation given by equation (3) yields a very ac- 
curate numerical representation for derivatives of /. 
Spectral methods therefore satisfy the requirements 
for approximating spatial derivatives in the equations 
governing a chemically reacting flow field. 

Two second-order finite-difference discretizations 
of the spatial derivatives are also developed, both 
to demonstrate the advantages offered by the higher 
order scheme and to provide benchmark results with 
more conventional approaches. In the first approach, 
second-order central finite differences are used to 
discretize the spatial derivatives. In the second 
approach, first-order forward and backward finite 
differences are used in combination with a predictor- 
corrector temporal discretization to yield a second- 
order method in space. 

With the approaches described above for tempo- 
rally and spatially discretizing the governing equa- 
tions, three numerical algorithms, one using spec- 
tral methods and two using finite-difference schemes, 
are developed for solving the equations governing 
a chemically reacting flow. The spectral algo- 
rithm (ref. 30) employs a two-stage partial implicit 
Runge-Kutta scheme for integrating the equations 
in time (ref. 23) and a Chebyshev spectral collo- 
cation method for computing spatial derivatives in 
the equations. The first finite- difference algorithm 
uses a partial implicit Adams-Moulton scheme to in- 
tegrate the equations in time and central finite dif- 
ferences to integrate the equations in space. The 
second finite- difference scheme employs a partial im- 
plicit MacCormack predictor-corrector scheme to in- 
tegrate the governing equations in time and space 
(refs. 23, 24, and 31). Computer programs have been 
written to apply these algorithms to the solution of 
reacting flow problems (ref. 30). The codes are lim- 
ited in this section to quas i- one-dimensional inviscid 
flows with hydrogen-air reaction, which is appropri- 
ate for the development and evaluation of the al- 
gorithms. Chemical reaction is represented in the 


programs with a finite-rate chemistry model, and a 
real gas thermodynamic model is employed. 

2.1 Governing Equations 

The quasi-one- dimensional Euler equations in 
conservation law form with multiple species under- 
going chemical reaction are (ref. 32) 


du dF __ 

~m + d^ + li = 0 

(4) 

where for i = I, 2, . . . , N s — 1, 


U = {pA y puA,pe 0 A,pf t A} T (i =1,2,.. 

...N,- 1) (5) 

F = [puA^pu 2 A + pA, puh 0 A y puf t A} T 

(6) 

H= 

(7) 

and 


fT 2 N * 

h 0 ~ 1 c p dT ^ (^r) l ft 

Jt R t=i 

(8) 


(9) 


where (H^)i is the reference enthalpy of species i 
at the reference temperature T# = 0 K (ref. 33). If 
there are N s chemical species, then i = 1,2,..., (N s — 
1) and (N s - 1) equations must be solved for the 
species / t . The final species mass fraction fjy s can 
then be found by conservation of mass since 

N s 

£/.- 1 

i = 1 


2.2 Chemistry Model 

The chemical reaction of hydrogen and oxygen is 
modeled here with the global finite-rate hydrogen- 
air chemistry model of Rogers and Chinitz (ref. 34). 
This model adequately represents the chemical reac- 
tion taking place in the problems to be considered 
in this chapter, and it also produces an extremely 
large disparity in the time scales in the problems. 
This phenomenon allows the ability of the numeri- 
cal algorithm to deal with resulting stiffness to be 
demonstrated. 

The Rogers-Chinitz model assumes that the over- 
all reaction of hydrogen and oxygen takes place 
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through two reactions, the first resulting in the for- 
mation of the hydroxyl radical, and the second com- 
bining the hydroxyl radical with hydrogen to form 
water. (More general models are needed, however, to 
properly represent the ignition stage of hydrogen-air 
reaction. These models are used in the more physi- 
cally detailed work described in section 3.) The re- 
actions are given by 

*/i 

(1) H 2 + 0 2 20H 

k bi 

k f 2 

(2) 20H + H 2 * 2H 2 0 

k b2 

where kf is the forward reaction rate and kf, is the 
reverse reaction rate. The reverse rate can be found 
given the forward rate and equilibrium constant K 
for each reaction, as 

k b = k } /K (10) 

The forward reaction rates are computed from the 
modified Arrhenius law, 

k ft = A l T N 'e~ E 'l R ° T (11) 

for each reaction i. For the Rogers-Chinitz model, 

the rates are given by (ref. 34) 


kf 1 

= A\T~ 10 e~ 4m5 / R ° T , 

cm 

mole-sec 

(12) 

k f2 

_ £ j <- 13 e -42500 /R°T 

cm 6 

(13) 

mole 2 — sec 


where 

A x = (8.9170 + 3L J 3 ? - 28.95) 10 47 
a 2 = ( 2.0 + - 0.8330) 10 64 

and 

K x = 26.164 e ~ 8992 / r 

K 2 = 2.682 x 1(T 6 T e 69415/T 

Knowing the reaction rates for the reactions defined 
by (1) and (2), the production of the four species 


present in the model can be found from the law of 
mass action. For a general reaction 

(i = 1,2, , Nji) 

>= 1 k b t J =1 

the law of mass action states that the rate of change 
of concentration of species j by reaction i is given by 
(ref. 35) 

&), = w- - <i) (*/. n c 1 is - fi 

V J=1 3 = 1 I 

(14) 

The rate change in concentration of species j by 
all Nft reactions is then found by summing the 
contributions from each reaction, 

< 15 > 

1=1 

Finally, the production rate of species j is found from 

Wj = CjMj (16) 

Applying the law of mass action to the global model, 
reactions (1) and (2), gives (ref. 34) 

Co 2 = -fc/iC H2 Co 2 + k bx C 2 on (17) 

^h 2 0 = 2 (fc /2 C& H <?H 2 - *W^h 2 o) (18) 

Ch 2 = Cq 2 - 2^H 2 0 (19) 

C OH = -(2Co 2 +Ch 2 o) (20) 

The source terms for the last i equations in equa- 
tion (4) can now be determined, as a function of the 
dependent variables, by application of equation (16). 

2.3 Thermodynamics Model 

The specific heat at constant pressure, c p , is 
nearly a linear function of temperature for each 
species present in the flow field (H 2 , 02, OH, H 2 O, 
N 2 ) over the range of temperature being considered 
in this section. To simplify the analysis, c p versus 
temperature data (ref. 33) for each species i is there- 
fore fit with 

c Pl (T) = a t T + hi (21) 
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where a and b are constants. A mixture specific heat, 
c p , can then be defined by weighting over the species 
i as 

N s 

Cp = JS ^2i C P x fi ( 22 ) 

i=l 

The total enthalpy of the mixture, made up of the 
five species, is given by 

H = £ U dT + («?)■)+ y (23) 

Putting equation (22) into (23) and integrating gives 


_ s ( n ,rp 2 

-5 "(t 


— + 6 i r + (flj) i + T 


Finally, the mixture gets constant i? is found by 
weighting the individual gas constants over the 
species i as 

Ns 

R = 'Z R ifi (25) 

Equations (22), (24), and (25) can then be used to 
define all other required thermodynamic variables. 

2.4 Chebyshev Spectral Method 

2.4.1 Spatial discretization. The Chebyshev spec- 
tral collocation method (ref. 28) is used to define the 
derivatives dF/dx in equation (4). To define dF/dx , 
F is expanded in terms of the Chebyshev polynomials 

T n (x) = cos(ncos~ 1 x) (26) 

in the truncated Chebyshev series 


F(x) = £ F n T n (x) 


where F n represents the expansion coefficients of the 
series. To form a range on x, the change of variables 


x = cos 6 


(0 < 0 < tt) (28) 


is introduced. Putting equation (28) into (26) and 
introducing the resulting expression into (27) gives 


F(x) — E F n cos(n0) 


a Fourier cosine series. To discretize equation (29), 
a set of collocation points Xj is defined by 


(j = 0,1, 2,... ,7V) (30) 


Xj = cos I — 


and the discrete form of equation (29) becomes 


Fj = F(xj) = ^2 F n i 


The inverse of equation (31) gives F n and can be 
found as follows. First, multiply equation (27) on 
both sides by T m (x) and weight (1 — x 2 ) 1 / 2 , and 
then integrate over the interval [-1,1]. This gives 

J' F(x) (l-x 2 )‘ 1/2 T m (x) dx 

= E F " T (i -x 2 )“ 1/2 r n (x)r m (x) dx 

n =0 J ~ 1 

Making the transformation x = cos 6 yields 
f' (l - x 2 )“ 1/2 T n (x)T m (x) dx 


rO 

= — I cos(n0) cos(m0) 
J 7T 


— 2 


where 


Therefore 


2 n = 0 or n = N 
1 1 < n < (N — 1) 


fn = ^L l ( 1 ~ x2 ) 1/2 F (x)T n (x)dx 


Again let x = cos 6 and 


^ 2 /' n7r 

F n = — — / F (6) cos (n0) dinO ) 

7rc n n Jo 

To generate a discrete set of Chebyshev coefficients, 
the trapezoidal rule of integration 

N i N 

I = ^2 ~ o JZ (9n + 9n+ 1) 
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is employed, where I is the value of the integral, h is 
the integration step size, and g is the integrand. The 
expression for discrete values of F n then becomes 


N- 1 


F n = 


7 vc n n 2 N 


[FjTn (^j) + F J + iT n (^j -hi)] 


3=0 


1 

c n N 


N — 1 N 

FjTn (xy) + y; FjTn (xj) 


3=0 


3 = 1 


1 2F 0 Tn(x 0 ) , o r / x , 2F N T n (x N ) 

■ 2 + 2 F > Tn <** J + 2 

J=1 


2 

CnW 


N 

'Z'ij'VjTniXj) 


3 = 0 


Putting equation (35) into (34) and algebraically 
manipulating the resulting expression gives 


N r ^( 1 ) N p(l) 

06) 


n=l 


2n 


Tl= 1 


Introducing equation (33) into (36) and simplifying 
then results in 


2nF„ = C n - 1 F^ 1 -F^ 1 


(37) 


an expression for given F n . The procedure for 
finding F^ is initialized by setting 


f?(l) n 

* N+l ~ 0 


Returning for consistency to the trigonometric form 
gives the discrete Chebyshev coefficients as 


N 


F n = 


Nc n 


Es lF V cos 

3=0 


/ nirj 

V n 


(32) 


Examination of equations (31) and (32) shows that 
F n can be efficiently evaluated using the fast Fourier 
transform (ref. 36). 

Next, F is differentiated in equation (31) with 
respect to x, giving 


N 

F'(x) = E F nT'(x) (33) 

n— 1 

A form of equation (33) without derivatives of the 
Chebyshev polynomials is preferred, so equation (33) 
is rewritten in terms of another series 

F'(z) = E Vn^Tnix) (34) 

n=0 


and then the coefficients of the two series are com- 
pared. The following recursion relation exists be- 
tween the Chebyshev polynomials and their deriva- 
tives (ref. 28) 


where 


rpf rpf 0 

1 n+l 1 Tl— 1 ^ rp 

r — r — TT'-in 

n+l n — 1 O fi 


2 n — 0 

1 n > 1 


(35) 



and then solving for F^ j through Fq 1 ^ by back sub- 
stitution (ref. 28). Then, knowing all F^\ the re- 
quired spatial derivatives of F can be calculated from 
equation (34). This procedure can again be done ef- 
ficiently with the fast Fourier transform (FFT). 

When the number of computational grid points 
to be used in a calculation is less than 60, it becomes 
more efficient to abandon the use of FFT’s and for- 
mulate an alternative method for spectrally comput- 
ing derivatives of F. The derivative is first written 
discretely as 


N 

F ; (x*:) — E D k jV(xj) (38) 

3=0 

where D^j is a matrix (termed the Chebyshev ma- 
trix) that must be found. An interpolant of F(x) 
at any point x must then be constructed. Following 
reference 37, the following polynomials are chosen 

, , (i-* 2 KM<-i > )+1 

0 > fx> - ljW)x - Xj) 

where Cj is defined the same as c n on page 9. The 
TVth degree interpolation of F(x) is then given by 

N 

F(x) = E 9j(x)F( Xj ) 

3=0 

To find F'(:r), the above expression is differentiated 
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where K™ is the Jacobian of H;, <9H/dU. Putting 
equation (41) into (40), simplifying the resulting 
equation, and then rewriting in delta form gives 

[/ + At K? } AU^ +1 = -At (^)" + H” 

' ' l sp 

(42) 

where [I] is the identity matrix and AU™ +1 = 

— U» Examination of equation (42) shows 
that the bracketed term on the left-hand side is a 
block-diagonal matrix, the blocks being n by n sub- 
matrices with n the number of equations in equa- 
tion (4). Since the matrix in equation (42) is diago- 
nal, equation (42) is the most easily solved for AU 
by inverting the blocks, i.e., 

AU” +1 = -At [/ + At K?}- 1 R" (43) 
where [ ] _1 represents a block invert, and 


The required derivatives of F can now be found by 
using relations (39) in equation (40). Because of the 
form of equation (38), this approach for computing 
F ; is often termed the direct matrix method. 

2.4.2 Temporal integration. Once values for dF jdx 
and H are determined as described above, there 
remains a system of ordinary differential equations in 
time that must be solved for the dependent variable 
vector U. The equations are integrated in time using 
a two-stage Runge-Kutta technique. The algorithm 
is developed as follows. 

Equation (4) is first discretized in time as noted 
above, giving 


R ?=(?■ 


+ H? 


is the steady-state residual vector. The two-stage 
Runge-Kutta technique is then applied to equa- 
tion (42), yielding the following predictor-corrector 
formulas: 


Predictor: 


AU" +1 = —At [I + At K ”}- 1 Rf 
U n+1 = U? + AU n+1 

t l l 


(44a) 


Uf +1 = U" — At (^)" + H" +1 + 0(Af) 2 

' ' *sp 

(40) 

where n is the old time level, n + 1 is the new time 
level, and sp indicates that the spatial derivatives are 
computed spectrally. Note that the source term is 
written implicitly as previously discussed to counter 
the potential effects of stiffness that may be encoun- 
tered in the governing equation system. The vector 
H n + 1 is then expanded in a Taylor series in time. 

H n+1 = H? + At (^ )" + O(At) 2 
or 

H” +1 = H” + K? (U t n+1 - U t n ) + O(At) 2 (41) 


Corrector: 


AU" +1 = —At [/ + At iff]" 1 R? +1 

. i (44b) 

U" +1 = Uf + - f AU n+1 + AU n+1 ') 
t 1 2 V 1 1 ) 

Starting with initial conditions for U, equations (44) 
are used to advance the solution from time level n 
to time level n - hi. The process is continued until 
steady-state conditions, defined as a reduction of 
10 orders of magnitude in the steady-state residuals, 
are reached. 

The magnitude of the time step in equations (44) 
is chosen based on the physical time scales present 
at any given time in the solution. The fluid-dynamic 
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time step A tj can be shown numerically to be limited 
by the Courant condition, 


At j — 


Ax 

u\+c 


(45) 


The chemical relaxation time for a species i is given 
by (ref. 38) 

tc = — (46) 

Changes in this relaxation time are then given by 


Ate 


Mpfi) 

w. 


(47) 


since w t remains nearly constant over a time step. 
For accuracy, it is required that the chemical time 
step be chosen such that no change in specific mass 
fraction (pf t ) greater than 0.0001 occurs over that 
time step. The computational time step At is then 
chosen to be the minimum over all grid points of the 
fluid and chemical time step, i.e., 


A£ = min(A2y, A£ c ) (48) 


2.5 Adams-Moulton Finite-Difference 
Scheme 


2.5.7 Spatial discretization . Central finite dif- 
ferences are chosen to define the spatial derivatives 
dF/dx in equation (4) for use with the Adams- 
Moulton time-stepping scheme. The spatial dis- 
cretization of F then becomes 



F ? + l 


-F r 


2Ax 


+ o{A X y 


(49) 


Note that the finite-difference representation of spa- 
tial derivatives is local in nature, whereas the spec- 
tral method of section 2.4 represents these derivatives 
globally. 


2.5.2 Temporal integration . Again knowing values 
for dF/dx and H, the resulting system of ordinary 
differential equations must be integrated in time. 
Equation (4) is discretized in time by using the 
Adams-Moulton scheme (ref. 39) to yield 


where a is the degree of implicitness. Proceeding 
as was done in section 2.4.2, F n+1 and H n+1 are 
expanded in a Taylor series in time to give 

F n+ 1 =F n + J n (u n+1 -U n ) +0(A<) 2 (51) 

H n+1 =H n + /f n (U n+1 -U") + 0{At) 2 (52) 

where J n is the Jacobian of F, <9F/dU. Putting 
equations (51) and (52) into (50), simplifying the 
resulting equation, and rewriting in delta form then 
gives 



where again AU n+1 = U n+1 — U n and fd indicates 
that the spatial derivatives are computed using cen- 
tral finite differences. In discrete form, equation (53) 
becomes 

,+ “ a, (^ +k ") 

(54) 

where 6 is a spatial central difference operator oper- 
ating on J and AU, and R is the steady-state resid- 
ual given by 


AU" +1 = -At R? 


F? 

R? - 1+1 


F n 

i—l i Tin 

2 Ax ^ H * 


(55) 


The bracketed term on the left-hand side of equa- 
tion (54) is a block tridiagonal matrix. This system 
can be solved using the Thomas algorithm (ref. 40). 
To apply that algorithm, equation (54) is rewritten 
as 

A n AU «+1 + B n AlJ n+l + Qn AlJ n+l = D” (56) 
where 

n _ a At 

“ 2Ai J ’- 1 

B? = {I + ocK™ At) 

n Q jn 

" 2~Ax ,+1 


U n+1 = U n - At 




n 

+ H n 


+ a 



+ H n+1 


+ 0{At) 2 (50) 


D" = -At R? 

It is then assumed that equation (56) can be written 
in upper triangular form as 

AU" +1 = E™ + F" AU^, 1 (57) 
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Putting equation (57) evaluated for AU" +1 into (56) 
and manipulating then gives 

AU," +1 = (B” + CfF^y 1 (d? - C?E? +l ) 

-(*? + C?F? +1 )~ l A? AU?+/ (58) 
Comparing equation (57) with (58) then yields 

E? = {B? + C?F t n +i y' (d? - C?E? +l ) (59) 

F? = -(B? + C?F 1 n +i y 1 A? (60) 

Once boundary conditions have been established 
(section 2.7), values for E™ and FJ 1 can be found 
by back substitution. Then, knowing these values, 
AU” +1 can be found by forward substitution from 
known values of AU^ 1 . Starting with initial con- 
ditions for U, equation (57) is used to advance the 
solution from time level n to n + 1. The process is 
continued until steady-state conditions are reached. 

The magnitude of the time step used to evaluate 
the coefficients in equation (57) is again chosen as 
in the spectral algorithm based on the physical time 
scales present in the problem. This choice is nec- 
essary to preserve the real-time accuracy of the so- 
lution. With the Adams-Moulton method, however, 
the time step chosen can be significantly larger than 
the time step based on physical time scales, since the 
method can be made fully implicit with proper choice 
of the implicitness factor a. The Adams-Moulton 
method is still attractive for real-time studies be- 
cause of its effective damping of high-frequency com- 
ponents present in the solution at early times. 

2.6 MacCormack Finite-Difference 
Scheme 

2.6.1 Spatial discretization. The MacCormack 
finite-difference method (ref. 31) is a predictor- 
corrector scheme of the Lax-Wendroff type. First- 
order forward differences 


the method becomes second-order accurate in space 
(ref. 31). It should be noted that the method can 
be made nearly symmetric by alternating the spatial 
differencing in the predictor and corrector steps with 
each succeeding time step, i.e., forward differences in 
the predictor and backward differences in the correc- 
tor on the first time step and backward differences 
in the predictor and forward differences in the cor- 
rector on the second time step, etc. The symmetric 
algorithm is applied in this work. 

2.6.2 Temporal integration. With a redefinition 
of the steady-state residual, the temporal integrator 
in the MacCormack scheme is identical to that em- 
ployed with the spectral spatial discretization, equa- 
tions (43) and (44), in Section 2.4.2. For the first 
time step, the predictor step residual is given by 


F? 

R n = i+l 


F? 


Ax 


+ H” 


(63) 


and the corrector step residual is given by 


p*n pn 

R” = (64) 

For the second time step, the residual definitions are 
alternated, and the process is continued until steady- 
state conditions are achieved. 

2.7 Initial and Boundary Conditions 

Governing equation (4) is hyperbolic and requires 
initial conditions at each point to start the calcula- 
tion and boundary conditions at the inflow boundary. 
Initial conditions are computed by first specifying an 
inflow Mach number and estimating an outflow Mach 
number. The interior Mach number distribution is 
then assumed to have a spatial variation that is lin- 
ear. The total pressure and total temperature are 
assumed to be constant throughout the domain. Fi- 
nally, the initial flow is assumed to be isentropic, so 
that isentropic relations can be used to compute the 
static pressure and temperature; these conditions are 
found from 



are used in the predictor step of the algorithm, and 
first-order backward differences 



are used in the corrector step. When these dif- 
ferences are summed in a predictor-corrector pass, 



(65) 

II 

<£!«■ 

(66) 


Knowing the static temperature, static pressure, and 
Mach number, the velocity distribution can be com- 
puted, and the density distribution can be found 
from the equation of state. Since the inflow bound- 
ary flow remains supersonic, boundary conditions are 
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specified there by holding conditions fixed at their 
initial values. 

The Adams-Moulton and MacCormack finite- 
difference schemes also require numerical boundary 
conditions at the outflow boundary. The spectral 
algorithm requires no such outflow boundary condi- 
tions since spatial derivatives can be defined at the 
outflow boundary in the same manner as is done at 
interior points. Outflow boundary conditions are de- 
fined for the finite-difference codes by using a second- 
order extrapolation formula. This formula is de- 
rived by writing a Taylor expansion to second or- 
der of the conserved variable vector U at the outflow 
boundary: 

U N = U ^Axjv-i+OtAx) 2 

In a spatially discrete form, this becomes 

Ujv = u N-l + UjV A 1 ~ + O(Ax) 2 

or 

U N = 2U N -i -Ujv-2 (67) 

since 

Ax N -l = AjjV-2 

Equation (67) is used directly in the MacCor- 
mack algorithm to define explicitly the numerical 
outflow boundary conditions. Boundary conditions 
are defined implicitly in the Adams-Moulton algo- 
rithm, however. To satisfy equation (67), equation 
(57) must be rewritten at the outflow node N to in- 
clude the N - 2 node. This is done by introducing a 
new coefficient G* y such that 

AU"/ 1 =E% + F% AU"t\ + Gn AU n N + _\ (68) 

and requiring that E jy = 0, Fjy = 2, and G ] y = 
— 1. This completes the definition of all required 
physical and numerical boundary conditions for the 
algorithms. 

2.8 Results 

Three numerical algorithms have now been devel- 
oped for solving the equations governing an inviscid 
chemically reacting flow; these algorithms were used 
to calculate the reacting flow in a rapid expansion 
supersonic diffuser. A rapid expansion diffuser was 
chosen such that high concentration gradients existed 
near the inflow boundary, providing a rigorous test of 
the methods. The comparison also allowed a demon- 
stration of performance of the high order accurate 


spectral method on grids that were quite coarse com- 
pared with the grids required for equivalent accuracy 
using the 2 finite-difference methods. The two finite- 
difference codes were also compared with each other 
to determine their relative accuracies and efficiencies 
when used to compute the test problem. 

The rapid expansion diffuser is shown in figure 3. 
The diffuser is 2 units long, has an initial cross- 
sectional area of 0.79 and a final cross-sectional area 
of 3.14. The diffuser wall is defined, as noted, by a 
shifted sinusoid. Flow is introduced to the diffuser 
at a Mach number of 1.4, a velocity of 1230 m/s, a 
temperature of 1900 K, and a pressure of 0.081 MPa. 
The chemical composition of the inflow is defined to 
be a three-tenths stoichiometric mixture of hydrogen 
fuel and air. 

Starting from the initial state described above, 
the governing equations are solved, using the three 
algorithms in a time consistent manner, until steady- 
state conditions are reached. A comparison of the 
spectral and finite-difference methods, as shown by 
the history of the chemical species, is given in fig- 
ures 4 through 6 for H 2 , O 2 , OH, and H 2 O, respec- 
tively. Results are presented at the first grid point 
interior to the inflow boundary, where the flow field 
and species gradients are a maximum. Agreement 
between the Runge-Kutta spectral code and the two 
finite-difference calculations is excellent in all cases. 

Next, spatial results from the methods are com- 
pared once steady-state conditions have been 
reached. The finite-difference solutions required 
101 grid points before a grid independent solution, 
defined as a graphically imperceptible difference in 
the steady-state result between the present grid and 
next coarser grid, was attained. Calculations using 
the Runge-Kutta spectral code were carried out on 
17- and 9-point grids. Steady-state results for the 
methods are given in figures 7 through 12. Figure 7 
shows the axial velocity distributions in the diffuser. 
The 17-point spectral solution and the 101-point 
finite-difference solutions agree quite well throughout 
the diffuser. The 9-point spectral solution slightly 
overpredicts the velocity near the inflow boundary, 
but agrees well throughout the remainder of the dif- 
fuser. The overprediction is likely to be due to the 
failure of the coarsest spectral grid to predict ade- 
quately the high gradients that exist at the beginning 
of the diffuser. Temperature distributions, given in 
figure 8, follow similar trends, with the 17-point spec- 
tral solution agreeing well with the difference calcu- 
lations, and the 9-point solution also agreeing well, 
except near the inflow boundary. Identical trends 
also occur when axial pressure distributions are 
compared in figure 9. 
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Axial species distributions computed by the 
methods are given in figures 10 through 12. Pre- 
diction of the H 2 mass fraction by the spectral 
method with 17 grid points agrees well with the 
finite-difference solution throughout the diffuser, as 
can be seen by examining figure 10. The 9-point spec- 
tral solution underpredicts the H 2 mass fraction near 
the inflow boundary, again due to the high spatial 
gradient in /// 2 there, but agreement again becomes 
good away from the inflow boundary. The spatial 
distribution of O 2 mass fraction is given in figure 11. 
The gradients are not as large for this species since 
O 2 is in excess, and both 17- and 9-point grids agree 
well with the finite-difference solution. The steady- 
state species distributions for OH and H 2 O are given 
in figure 12. The spatial gradients are again high for 
both species near the inflow boundary, and trends 
similar to those for H 2 are repeated here. Agree- 
ment is again quite good when comparing the 17- 
point spectral and finite- difference results. The 9- 
point spectral solution still underpredicts gradients 
near the inflow boundary, however. 

A final comparison of methods can be made in fig- 
ure 13, which shows the rate of reduction of steady- 
state residual with iteration count for each algo- 
rithm at the first interior grid point. Since the 
17-point Runge-Kutta spectral and the 101-point 
finite-difference calculations yield comparable accu- 
racy and have the same minimum spatial step size, it 
is reasonable to assess the relative efficiency of the 
methods by using the results given in this figure. 
Note that the residual reduction rate by the spec- 
tral code is significantly greater than that provided 
by the finite-difference codes. The maximum residual 
(at any grid point) is reduced with the spectral code 
by 10 orders of magnitude in only 2400 iterations. 
The Adams-Moulton finite-difference code requires 
4000 iterations to achieve the same level of residual 
reduction. The MacCormack finite-difference code 
is only able to achieve a 2- to 3-order-of-magnitude 
reduction in steady-state residual because of its in- 
consistent residual definition between predictor and 
correction steps. (Recall that forward spatial differ- 
ences are used in the predictor, and they are alter- 
nated with backward differences in the corrector.) 
Even with this deficiency, however, the MacCormack 
method is able to achieve an acceptable level of ac- 
curacy in comparison with the Adams-Moulton and 
spectral schemes, as can be seen from the previous 
results. 

To achieve a fair comparison of the three algo- 
rithms, the convergence history discussed above must 
be combined with the computational grid needed to 
achieve the required accuracy and the computational 
time required per time step for each scheme. The 


Runge-Kutta spectral code on the 17-point grid re- 
quired 644 CPU seconds to meet the established con- 
vergence requirement. The Adams-Moulton code on 
the 101-point grid required 1706 CPU seconds to also 
meet the convergence requirement. As noted before, 
the MacCormack code did not meet the convergence 
criteria, but it did achieve an acceptable level of accu- 
racy for a steady-state solution on the 101-point grid 
after 4000 iterations. The code required 876 CPU 
seconds to reach steady-state conditions. 

Based on the above results, algorithms to be ex- 
tended to two-dimensional flows were chosen. The 
Runge-Kutta spectral method was an obvious choice 
because of both its high accuracy and its excellent 
computational efficiency. The MacCormack method 
was computationally more efficient than the Adams- 
Moulton scheme, but it was unable to achieve as 
high a degree of steady-state residual reduction. One 
other fact, not apparent in the previous calculations, 
must be considered in this comparison, however. The 
Adams-Moulton scheme resulted in a system of equa- 
tions that contained block tridiagonal structure. The 
MacCormack scheme resulted in a system of equa- 
tions that contained only block diagonal structure if 
the system was stiff, and no left-hand-side matrix at 
all if the system was not stiff. The work required 
to solve a block tridiagonal system varied with 3/V 3 , 
where N is the number of equations. The work nec- 
essary to solve a block diagonal system increased 
with jV 3 , and the work to solve the system with- 
out a left-hand-side matrix increased with N. It was 
found in section 3 that when detailed (as opposed to 
global) chemistry systems were used to model super- 
sonic reacting flows, the resulting system of equations 
was not temporally stiff. When the points described 
above were considered in this light, the MacCormack 
algorithm became the preferred finite-difference al- 
gorithm of those considered for extension to two- 
dimensional flows. 

3. Multidimensional Chemically 
Reacting Flows 

In the previous chapter, three algorithms were de- 
veloped for the study of inviscid quasi-one- 
dimensional, supersonic, chemically reacting flow. 
From those algorithms, two were chosen for exten- 
sion to two-dimensional, viscous, supersonic, chemi- 
cally reacting flow. Those extensions are carried out 
in this section. Additionally, considerably more de- 
tailed chemistry and thermodynamics models are de- 
veloped here for the programs. Finally, to include the 
effects of diffusion of momentum, energy, and mass, 
kinetic-theory-based diffusive transport models are 
developed and incorporated into the programs. De- 
tails of these models are given in the following 
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section. They are discussed following a statement 
of the general system of equations governing two- 
dimensional, viscous, chemically reacting flows. 

3.1 Governing Equations 

The two-dimensional, Navier-Stokes, energy, and 
species continuity equations governing multiple 
species undergoing chemical reaction are given by 
(ref. 35) 
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The diffusion velocities are found by solving (ref. 35) 
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Note that if there are N s chemical species, then 
i — 1 , 2 ,..., (7V S — 1) and (N s — 1) equations must 
be solved for the species / t . The final species mass 
fraction //y can then be found by conservation of 
mass since 

N s 

£/.-i 

i= 1 


3.2 Thermodynamics Model 

To calculate the required thermodynamic quanti- 
ties, the specific heat for each species is first defined 
by a fourth-order polynomial in temperature: 

^ = A { + B { T + CiT 2 + DiT 3 + EiT 4 (74) 
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The coefficients are found by a curve fit of the data 
tabulated in reference 33. Knowing the specific heat 
of each species, the enthalpy of each species can then 
be found from equation (71), and the total internal 
energy is computed from equation (70). 

To determine the equilibrium constant (required 
in section 3.3) for each chemical reaction being con- 
sidered, the Gibbs energy of each species must first 
be found. For a constant pressure process, c p /Rx 
from equation (74) is first integrated over temper- 
ature to define the entropy of the species, and the 
resulting expression is integrated again over temper- 
ature to obtain the following fifth- order polynomial 
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in temperature for the Gibbs energy of each species: 

^ = Ai(T -TlnT) - ^T 2 - ^T 3 - 0-T 4 - 0-T 5 + F t -G.T 
R 2 6 12 20 1 ’ 

(75) 

Coefficients F± and G t are again defined in refer- 
ence 33. The Gibbs energy of reaction can then be 
calculated as the difference between the Gibbs energy 
of product species and reactant species. 

AG R = E n i A 9i - n i ( 76 ) 

i=products i=reactants 

The equilibrium constant for each reaction is then 
found from (ref. 41) 


the law of mass action states that the rate of change 
of concentration of species i by reaction j is given by 
(ref. 35) 




(80) 


All third-body efficiencies are assumed to be equal 
to 1.0. The net rate of change in concentration of 
species i by reaction j is then found by summing the 
contributions from each reaction, 


Nr 

Q=£(c.), (81) 

3 = 1 



where An is the change in the number of moles when 
going from reactants to products. 

3.3 Chemistry Model 


Finally, the production of species i can be found from 

w { = CiM { (82) 

The source terms for the last i equations in (69) 
are now determined as a function of the dependent 
variables. 


In the present work, the finite-rate chemical reac- 
tion of gaseous hydrogen fuel and air is of concern. 
That reaction is modeled by a 9-species, 18-reaction 
model described in table I (ref. 42). Eight of the 
chemical species (H 2 , 0 2 , H 2 0, OH, H, O, H0 2 , 
H 2 0 2 ) are active, and the ninth (N 2 ) is assumed in- 
ert. The forward rate of each reaction j is given by 
the modified Arrhenius law 


3.4 Diffusion Models 

Models for the coefficients governing the diffusion 
of momentum, energy, and mass are now determined. 
The individual species viscosities are computed from 
Sutherland’s law, 


JL = /r\ 3/2 r 0 + s 

Ho \ T a ) T + S 


(83) 


k fj= A 3 T J N ( 78 ) 

Values for A, JV, and E are also given in table I. 
Knowing the forward rate, and using the equilibrium 
constant determined in the previous section, the 
backward rate can be defined by 

k bj = kf./Kj (79) 

Once the forward and reverse reaction rates have 
been determined, the production rates of the eight 
species can be found from the law of mass action. 
For the general chemical reaction 


where fi 0 and T 0 are reference values and 5 is the 
Sutherland constant. All three values are tabulated 
for the species in references 43 and 44. Once the 
viscosity of each species has been determined, the 
mixture viscosity is determined from Wilke’s law 
(ref. 45), 

= E % ( 84 ) 
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The species thermal conductivities are also computed 
from Sutherland’s law, 


k /T\ 3/2t o + S' 

k~o~\n) T + S' 


( 86 ) 


but with different values of the reference values k Q 
and and the Sutherland’s constant S'. These 
values are also taken from references 43 and 44. The 
mixture thermal conductivity is computed by using 
conductivity values for the individual species and 
Wassiljewa’s formula (ref. 46), 


_ t n 

i=1 1 + X- £ X A, 

J = 1 


suggested in reference 35, that the thermal diffu- 
sion coefficient D t is negligible compared with the 
binary diffusion coefficient. The solution of equa- 
tion (73) requires solving a simultaneous equation 
system, with the number of equations equivalent to 
the number of species present for each component of 
the diffusion velocity. It should be noted that for i 
species, however, the system of i equations defined by 
equation (5) is not linearly independent. One of the 
equations must be replaced by the constraint 

N s 

£p/;V t = 0 (92) 

1=1 

to make the system linearly independent. The re- 
sulting system of equations is solved for the diffusion 
velocities by using the Householder method (refs. 47 
and 48). 


where (\>- — 1.065<^y, and <p l3 is taken from 

equation (85). 

For dilute gases, Chapman and Cowling used 
kinetic theory to derive the following expression for 
the binary diffusion coefficient D tJ between species i 
and j (ref. 43): 

0.001858T 3 / 2 [(M,- + Mj) 1/2 

D i:j = 27 ; (° 8 ) 

P^ij^D 

Here, the diffusion collision integral fl/j is approxi- 
mated by 


n D = r— 0145 + (r* + o.5)‘ 


where 


T* = T/T e , 


Values of the effective temperature T e and effec- 
tive collision diameter a are taken to be averages of 
the separate molecular properties of each species, giv- 
ing (ref. 43) 

= j) ( 90 > 
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3.5 Solution of the Governing Equations 

Once the thermodynamic properties, diffusion co- 
efficients, and chemical production rates have been 
defined, the governing equations can be solved nu- 
merically. The finite-difference solution procedure is 
discussed in the next section, 3.5.1, and this discus- 
sion is followed by the development of the spectral 
solution scheme, described in section 3.5.2. 

3.5.1 Finite-difference solution method. To 

solve the governing equations (69) with the finite- 
difference scheme, the equations must first be trans- 
formed from the physical domain (x, y) in which they 
are written to an appropriate computational domain 
(£, r]). The equations are solved on a coordinate grid 
that is highly compressed in both x and y in the phys- 
ical domain near regions where high gradients exist. 
The grid is required to be uniform, however, in the 
computational domain to most readily maintain a re- 
quired order of accuracy. 

To transform the governing equations from the 
physical to the computational domain, fluxes F and 
G are first written in functional form and differenti- 
ated with respect to the computational coordinates. 
Therefore, given F = F (x,y) and G = G(x,y), and 
proceeding first with F, 

Ff = FjrX^ + Fy?/£ (93) 


Once the binary diffusion coefficients for all 
species combinations are known, the diffusion veloc- 
ities of each species can be computed from equa- 
tion (73). The diffusion velocity of each species is 
the species velocity due to all diffusion processes al- 
gebraically added to the convection velocity. When 
computing the diffusion velocities, it is assumed as 


F 77 — F x’X'T) H" F yVr) (94) 

Then, substituting F^ from equation (94) into equa- 
tion (93) and simplifying gives 

F, = ^ j (95) 
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where 


J = x^yrj - y^Xr, (96) 

is the Jacobian of the transformation. Proceeding in 
like manner for G gives 


_ G,,X£ - G^x,, 
Gy- J 


(97) 


layer. The compression function in the streamwise 
direction is given by 

£ = £o 

0 < £ < 1 
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Finally, substituting equations (95) and (97) into 
equations (69) gives the governing equations in the 
computational domain 


where 


<9U dF dG ^ 
~dt + df + Tfy “ 


(98) 


U = JU, H = JH 


F - yrj F - x v G 
G = X£G -y^F 

Here (x^, x^, y^, y^) are the transformation metrics 
that form the inverse Jacobian matrix, and J is the 
Jacobian of the transformation. The metrics can be 
computed numerically once the physical coordinate 
grid has been prescribed. 

To resolve large flow field and concentration gra- 
dients, the physical coordinate grid must be chosen 
sufficiently fine in those regions. For the mixing layer 
problems to be studied in this work, the grid must be 
highly refined in a direction transversely across the 
layer. Large streamwise gradients may also occur 
with movement along the layer, and grid refinement 
must also be allowed at those locations. The com- 
pression function of Thomas et al. (ref. 49) can be 
used to satisfy the refinement requirements in both 
the transverse direction and the streamwise direc- 
tion. The compression function in the transverse di- 
rection is given by 
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The degree of transverse compression is determined 
by f3 y , and rj 0 is the value of rj at which maximum 
compression occurs, i.e., the center of the mixing 


0 <£< 1 


where 



1 + (eP* - l) 

1 + (e~P* - 1) £ 0 


The degree of streamwise compression is determined 
by r, and £ 0 is the value of £ at which maximum 
compression occurs. 

Having now determined the nondimensional phys- 
ical domain (£, fj) from the computational domain 
(£>*7) by using equations (99) and (100), (£,r?) is 
then mapped onto the physical domain (x, y) by us- 
ing the algebraic-numerical coordinate transforma- 
tion of Smith and Weigel (ref. 50). This transforma- 
tion is given by 


*=x 2 (mv)+x 1 {i)[ i-ifo)] (ioi) 
y = *2(0 v{v) + *i( 0 [l - *?(*?)] (102) 

where {Xi>Yi) are the boundary points at y — 0, 
and (^ 2 , 12 ) are the boundary points at y = y max . 
The generality afforded by equation (101) that allows 
transverse coordinate lines to be skewed is not needed 
in this work. Therefore, X\ is chosen equal to X 2 and 
equation (101) simplifies to 


x — ^max (103) 

where £ is found from equation (100). The trans- 
formation metrics (x^, x^, y^, y^) are then found by 
directly differentiating equations (102) and (103). In- 
verse metrics (£ Xj £ y , rj x , rj y ) required for differentiat- 
ing terms within the flux vectors are then found by 
inverting the inverse Jacobian matrix, i.e., 


Vx 
Vy 

to form the Jacobian matrix of the transformation. 

It is sometimes advantageous to allow refinement 
of the physical grid in a point-by-point fashion. That 
option can be quite valuable for defining the stream- 
wise grid in the present work, and so such an option 
is provided by way of a simple modification of equa- 
tions (102) and (103). Rather than defining X\ and 
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X 2 with equation (100), the boundary points are de- 
fined manually in point-by-point fashion. Care must 
be taken that changes between successive points, Ax, 
are not too great or discretization errors can be in- 
troduced. Having defined the distribution of X\ and 
X 2 , the corresponding values of Y\ and Y 2 , and ff 
from equation (99), the required transformation met- 
rics can be found by numerically differentiating equa- 
tions (102) and (103) with respect to £ and rj . Once 
the choice for streamwise grid definition is made, 
all quantities required to describe the physical and 
computational domains are defined by equations (99) 
through (104). 

As noted in section 2, the governing equation sys- 
tem (69) can be stiff because of the kinetic source 
terms contained in the vector H and because of dif- 
fusive terms in the vectors F and G. Only the kinetic 
terms introduce stiffness in this work; spatial stiff- 
ness is controlled by the choice of grid. To deal with 
the stiff system, the approach used in references 23 
through 26 is again followed, and the kinetic source 
terms are computed implicitly. In a temporally dis- 
crete form, equation (69) then becomes 

n 

-H n+1 
(105) 

Following the approach used in section 2.4.2, H is 
linearized by expanding it in a Taylor series in time. 
Introducing this expression into equation (105), sim- 
plifying, and rewriting in delta form then gives 


modified MacCormack technique then becomes 
[/ - At K§] AuJ* 7 = - At r” (108) 

= U" + AU^ 

[/ - At K" +1 ] AU” +1 = -At R n+ ' (109) 

U” +1 = U n +0.5 (AU n+1 +AU n+1 ) 

where R represents a forward spatial difference of R 

and R a backward spatial difference. Stress terms 
are differenced in the conventional manner (ref. 31). 
Equations (108) and (109) are used to advance the 
solution from time n to time n + 1 and this process 
is continued until a desired integration time has been 
reached. 

The magnitude of the time step in equations (108) 
and (109) is again chosen based on the physical time 
scales present at any given time in the solution. 
These scales are defined in section 2.4.2; they are 
repeated here for convenience. The fluid-dynamic 
time step A tj can be shown to be limited by the 
Courant or viscous stability limit of the governing 
equation (ref. 31). The chemical relaxation time for 
a species i is given by (ref. 38) 



Changes in this relaxation time are then given by 



[/ - At K n ] AU n+1 = -At R n (106) 

where 



is the steady-state residual, I is the identity ma- 
trix, K n is the Jacobian of H with respect to U, 
(dH/dU), and AU n+1 = U n+1 - U n . 

Once the temporal discretization used to con- 
struct equation (106) has been performed, the re- 
sulting system is spatially differenced by following 
the approach of section 2.6 and using the unsplit 
MacCormack predictor-corrector scheme (ref. 31). 
This results in a spatially and temporally discrete, 
simultaneous system of equations at each grid point 
(refs. 23 and 24). Each simultaneous system is solved 
with the earlier noted Householder technique in com- 
bination with the MacCormack technique, which is 
then used to advance the equations in time. The 
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since W{ remains nearly constant over a time step. 
For accuracy, it is required that the chemical time 
step be chosen such that no change in mass fraction 
greater than 0.01 occurs over that time step. The 
computational time step A t is then chosen to be 
the minimum over all grid points of the fluid and 
chemical time steps. 


3.5.2 Hybrid Chebyshev spectral solution method. 

A hybrid Chebyshev spectral method has also been 
used to solve the governing equations (69) in this 
work. The spectral method, as discussed in section 2, 
is attractive for studies of reacting mixing layers be- 
cause it yields high numerical accuracy on relatively 
coarse grids. A highly accurate method is necessary 
for proper resolution of the large transverse gradi- 
ents that exist across the mixing layer. Gradients are 
not as large, however, in the streamwise direction of 
the mixing layer. A lower order method appeared 
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adequate in that direction. With these require- 
ments in mind, it was decided that transverse deriva- 
tives across the mixing layer should be computed 
spectrally, whereas finite differences were deemed 
appropriate in the streamwise direction. 

Spatial derivatives in the transverse direction 
were computed spectrally by using the direct 
Chebyshev matrix method developed in section 2.4.1. 
Required derivatives of the flux vector were com- 
puted at each grid point, given the distribution of 
the function Gy along the complete column of points 
which included point fc, i.e. , 
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employed. That function is given by 


V = 


1 — r; 2 


( 111 ) 


which maps the Chebyshev computational domain 
[-1,1] onto the physical domain [—00,00]. Maximum 
grid refinement occurs at y = 0, and the grid is cho- 
sen so that the mixing layer lies near this coordi- 
nate location. The function (3 y determines the de- 
gree of grid refinement. Equations (103) and (111) 
therefore complete the transformation from the com- 
putational to the physical domain. Elements of the 
inverse Jacobian matrix are again determined by di- 
rectly differentiating equations (103) and (111), and 
the Jacobian matrix is found by using equation (104). 

Having now defined the transformation required 
for the spectral method, the spatial derivatives are 
discretized as described earlier in this section. Once 
the spatial terms are differenced, there again remains 
a system of ordinary differential equations at each 
grid point to be integrated in time. Once the steady- 
state residual R? is redefined to reflect the change 
to spectral transverse derivatives, the procedure for 
temporally integrating the equations is identical to 
that carried out in equations (105) through (109). 
Introducing the new residual definition into equa- 
tions (108) and (109), the hybrid spectral algorithm 
is given by 


Streamwise spatial derivatives were again computed 
with the MacCormack finite-difference technique as 
applied in the previous section, 3.5.1. 

The governing equations must again be trans- 
formed from the physical (x, y) to the computational 
(£,77) domain, and the procedure described in equa- 
tions (93) through (98) is again employed. The 
streamwise compression function, equation (100), 
and the streamwise transformation, equation (103), 
are still equally valid, and they are retained. The 
streamwise grid can optionally be obtained in point- 
by-point fashion as before. The transformation in 
the transverse (spectral) direction must still be ca- 
pable of refining the grid at the center of the mixing 
layer. A different transformation is used, however, to 
preserve spectral accuracy when forming the trans- 
verse derivatives. Boyd (ref. 51) found that expo- 
nential mappings, such as the mapping employed in 
equation (99), gave poor performance. Calculations 
with alternate mapping functions indicated that, in 
general, the mapping function should decay more 
slowly than a function best describing the solution 
being sought. Recognizing the general form of the 
resulting mixing layer solution, an algebraic mapping 
function suggested by Boyd (ref. 51) was chosen and 
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where sp indicates that the derivative is to be evalu- 
ated spectrally. The time step At is again chosen by 
following the procedure described in section 3.5.1. 


3.5.3 Boundary and initial conditions. The govern- 
ing equations (69) require boundary conditions along 
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all four boundaries. For the problems to be consid- 
ered, the inflow boundary is always supersonic, so 
the velocities, static temperature and pressure, and 
species are specified and fixed there. The upper and 
lower boundaries always lie in the free stream, and 
therefore either the normal gradient of the preced- 
ing variables is required to vanish or the free-stream 
conditions are enforced along those boundaries. The 
gradient conditions not only satisfy the free-stream 
conditions, but also provide nonreflective conditions 
that pass disturbances through the boundary rather 
than reflect them back to the domain. The outflow 
boundary is also supersonic, and values of the veloci- 
ties, static temperature and pressure, and species are 
determined by extrapolation from upstream values. 
Finally, no slip boundary conditions (u = 0, v = 0) 
are used to specify velocity components along solid 
surfaces that occur in the physical domain. Addi- 
tionally, along these solid boundaries, adiabatic con- 
ditions ( dT /dy = 0) are assumed, the boundary- 
layer assumption on pressure ( dp/dy = 0) is cho- 
sen, and the walls are assumed to be noncatalytic 

( dfi/dy = o). 

The governing equations (69) also require a set 
of initial conditions. The equations are initialized 
by setting values of the velocities, static temperature 
and pressure, and species throughout the domain to 
the values chosen initially for boundary conditions 
at the inflow boundary. Having specified all required 
initial and boundary data, the equation is marched 
in time from the initial time to some final specified 
integration time. 

A general model for chemically reacting flow has 
now been developed, and the resulting governing 
equations have been defined in this section. Further, 
two numerical methods for solving these governing 
equations have been developed. In the following 
section, the governing equations are solved for several 
supersonic chemically reacting mixing layer cases, 
and the results are then discussed in light of the 
observations given in section 1 for such flows. 

4. Simulations of Reacting Mixing 
Layers 

By using the theory and solution procedure devel- 
oped in section 3, the chemically reacting flow field in 
a non-premixed laminar, supersonic, spatially devel- 
oping mixing layer is numerically simulated in this 
section. Two basic mixing layer cases are consid- 
ered. The first of those cases involves a hydrogen-air 
mixing layer with fuel and oxidant initially separated 
by a finite-thickness splitter plate. The second case 
considers a hydrogen-air mixing layer that has just 
begun to develop downstream of a splitter plate. The 


plate is not included in this calculation; the effects 
of the plate on the flow are retained, however. The 
first case is computed with only the finite-difference 
algorithm, and the results from that analysis are dis- 
cussed in the following section. The second case is 
computed with the hybrid spectral algorithm. Re- 
sults from that analysis are discussed in section 4.2. 

4.1 Simulations Using the Finite- 
Difference Algorithm 

The finite-difference algorithm has been applied 
to a non-premixed, spatially developing, laminar, su- 
personic, chemically reacting mixing layer. The con- 
figuration that is considered is described schemat- 
ically in figure 14, The overall domain is 5 cm 
high and 5 cm long. The height chosen places the 
boundaries well into the free stream, and the length 
allows initial development of the mixing layer. Ini- 
tially, hydrogen fuel and air are separated by a 0.5- 
cm-long splitter plate that is 0.02 cm thick and cen- 
tered at y — 2.5 cm. Downstream of the plate, the 
fuel and air mix and ignition occurs at some fur- 
ther distance downstream of the plate base. From 
that point, chemical reaction between the fuel and air 
takes place. For the problem being considered, cold 
gaseous hydrogen is introduced above the plate at 
Mach 1.5, a Reynolds number of 3700 based on plate 
thickness, a temperature of 293 K, and a pressure 
of 0.101 MPa (1 atm). Hot air is introduced below 
the plate at Mach 1.5, a Reynolds number of 731, a 
temperature of 2000 K, and a pressure of 0.101 MPa. 
These conditions result in an initial hydrogen veloc- 
ity of 1953 m/s and an air velocity of 1297 m/s; this 
yields a hydrogen- to- air velocity ratio of 1.5. 

By using the configuration and conditions de- 
scribed above, the mixing- layer flow field is marched 
in time from the specified initial conditions to the 
conditions existing at 0.1 ms. The solution is ob- 
tained on a spatial grid with 219 nodes in the stream- 
wise direction and 51 nodes in the transverse direc- 
tion. The grid is compressed in x near the trailing 
edge of the plate and highly compressed in y in the 
region of the mixing layer. The resulting flow field 
is described in figures 15 through 35, which give pic- 
tures of the flow at an instant in time. Figure 15 
shows a velocity vector field plot of the flow close to 
and on either side of the splitter plate and the de- 
veloping mixing layer downstream of the base of the 
splitter plate. (Velocity vectors are shown for only 
every four streamwise and transverse grid points in 
this region.) Expansions of both streams through 
Prandtl-Meyer fans can be observed at the trailing 
edge of the plate. The higher velocity hydrogen 
stream and the lower velocity airstream are appar- 
ent as is the wake flow downstream of the plate. The 
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development of the mixing layer with streamwise dis- 
tance can also be seen. Two regions of instability 
are also apparent in figure 15. The first region lies 
just downstream of the splitter plate approximately 
1.0 cm beyond the initial station. The second region 
lies well downstream at approximately 4.0 cm from 
inflow. The flow is relatively quiescent between these 
two regions. 

These instabilities can also be observed in fig- 
ure 16(a), which shows a plot of streamwise veloc- 
ity versus streamwise coordinate at several constant 
transverse stations that are well within the mixing 
layer. The oscillations are present along all three 
lines of constant y and are quite close in phase and 
magnitude, indicating that the instability is present 
in similar fashion across the layer. The oscillatory be- 
havior of the layer is quite typical of that seen numer- 
ically and experimentally in nonreacting layers and, 
at least in part, appears to be produced by a Kelvin- 
Helmholtz instability. To determine whether the in- 
stabilities and their locations were functions of time, 
several other times in the evolution of the layer are 
examined. These results are given in figures 16(b), 
16(c), and 16(d) for times of 0.09, 0.06, and 0.02 ms, 
respectively, beyond initiation of the flow. The in- 
stability near the splitter plate is present at all times 
that are given in the figures. The size of the region 
and the amplitude of the instability do not change 
with time. The location of the waves does change 
with time, however. At the latest time plotted the 
waves propagate downstream with increasing time, 
amplifying between x = 0.6 and 1.2 cm and damping 
beyond that streamwise station. The disturbance has 
essentially dissipated at 0.1 ms beyond 2.1 cm. At 
earlier times, however, the upstream instability prop- 
agates further downstream, reaching as can be seen 
in figure 16(d) the initiation of the second instability. 
With increasing time, though, the initial disturbance 
damps, and the central region of the flow between 
x = 2.1 and 2.8 cm becomes relatively quiescent. 

The initial velocity increase at x = 0.5 cm is also 
present at all times and is due to the expansion of 
hydrogen and air off the trailing edge of the split- 
ter plate. The velocity decrease that follows results 
from each gas being compressed by a recompression 
shock that turns the fluids to a nearly streamwise 
direction. In this region, in the wake of the plate 
just downstream of the splitter plate, the flow also 
separates. A recirculation bubble then forms and re- 
mains throughout the calculation. (The presence of 
the recirculating region can be seen in fig. 17; that is 
discussed later.) This separation is not stable with 
time; rather, it changes shape and position slightly 
with increasing time and acts as a destabilizing mech- 
anism for the flow in the wake downstream of the sep- 


aration. Changes in the position of the separation 
also cause the recompression shock to change posi- 
tion with time. The oscillatory motion with time 
of both the separation bubble and the recompres- 
sion shock thus appears to be the genesis of the first 
instability. It should be noted, however, that al- 
though the separation and recompression shocks ap- 
pear to be the tripping mechanism for the first insta- 
bility, the numerical method being employed suffers 
to some degree from Gibbs oscillations in the imme- 
diate neighborhood of the shock. These numerical os- 
cillations may also contribute to the initiation of the 
instability. 

The second instability is also present at all 
times shown in figures 16(a) through (d). Initially 
(fig. 16(d)), this region is influenced by the upstream 
disturbance. With increasing time, however, the re- 
gion preceding this instability becomes stable as can 
be seen by viewing figures 16(c), (b), and (a). It 
appears that the second region of instability repre- 
sents the onset of transition in the mixing layer. The 
amplitude of the disturbance grows with increasing 
distance downstream from the 2.8-cm station. There 
is also some growth in amplitude with time at any 
given streamwise station within the region of the 
instability. 

To examine the contribution to instability from 
heat release due to chemical reaction, the identical 
flow field is computed without reaction. Those re- 
sults are given for two times, 0.1 ms and 0.02 ms, 
in figures 16(e) and (f), respectively. By comparing 
figures 16(e) and (f) with 16(a) and (d), it can be 
seen that both instabilities still remain without heat 
release. The upstream disturbance, in fact, appears 
essentially unaffected by reaction. The effect of heat 
release on the downstream disturbance also appears 
quite mild at early times (figs. 16(d) and (f)), but 
there is a marked effect at 0.1 ms, as can be seen by 
comparing figures 16(a) and (e). The amplitude of 
the disturbance is consistently greater without chem- 
ical reaction. This result is consistent with the find- 
ings of references 15 and 20 for subsonic flow, which 
showed mixing is retarded by heat addition. 

Another view of the streamwise development of 
the velocity field is given in figure 17, which shows 
u profiles as a function of the y coordinate at four 
(x = 0.51, 1.0, 3.0, and 5.0 cm) streamwise stations. 
The initial profile shows that there is a recirculation 
region with negative streamwise velocities near the 
trailing edge of the splitter plate. A velocity defect 
in the wake continues to exist downstream at the 
1.0- and 3.0-cm stations. A developed mixing layer 
profile is apparent once the 5.0-cm station is reached. 
An overall view of the streamwise velocity is given 
by using a contour plot in figure 18. The regions 
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of instability and the development of the mixing 
layer are consistent with figures 16 and 17, but 
they can be viewed in a more realistic sense when 
shown in two dimensions. The velocity contours 
are compared with two-dimensional contours of other 
primitive variables later in this discussion. 

A plot of static temperature versus streamwise 
coordinate for several constant transverse stations is 
given in figure 19. The y coordinates are identical to 
those given in figure 16. The instabilities present in 
the velocity plots of figure 16 are consistent in loca- 
tion with those of the temperature field. The ampli- 
tudes of the oscillations in temperature are greater, 
however, because of a significantly greater temper- 
ature difference between hydrogen fuel (293 K) and 
air (2000 K) as compared with the velocity difference 
between the fuel stream and airstream. The first 
disturbance is significantly more pronounced along 
the y — 2.5 cm coordinate line where cold fuel and 
hot air are initially in contact as compared with the 
y — 2.46 cm and y — 2.54 cm lines where no mix- 
ing occurs. The second disturbance is markedly more 
pronounced along the y = 2.46 cm and y = 2.5 cm co- 
ordinate lines as compared with the y = 2.54 cm line, 
indicating that thermal mixing is occurring mainly 
below the location of the splitter plate at y = 2.5 cm. 

A plot of temperature profiles versus the trans- 
verse coordinate at four streamwise stations (x = 
0.51, 1.0, 3.0, and 5.0 cm) is given in figure 20. These 
stations are the same as those used in figure 17. The 
development of the temperature profile with increas- 
ing streamwise distance can be seen in the figure. 
Initially (x — 0.51 cm) there is some cooling in the 
base region beyond the plate because of expansion 
of the fuel and air off the plate and because of the 
endothermic reactions associated with ignition early 
in the development of the layer. Well downstream 
at x — 5.0 cm, the temperature profile is well devel- 
oped, and there are temperature increases on either 
edge of the layer associated with the exothermic re- 
actions that are taking place. 

By comparing figures 20 and 17, it can be seen 
that the temperature profiles at each x-station are 
consistently broader than the streamwise velocity 
profiles. This is also apparent by comparing the plot 
of temperature contours in figure 21 with the veloc- 
ity contours given in figure 18. This result is con- 
sistent with the discussion and experimental obser- 
vations described earlier from references 1 through 
20, and in particular reference 16. Vortical struc- 
tures are present in the mixing layer, and the ex- 
istence and growth of these vortices influence the 
growth and reaction in the mixing layer. The vor- 
tical character can be seen in figure 22, which gives 
the vorticity distribution in the mixing layer. Chemi- 


cal reaction takes place not only in the interior of the 
mixing layer, but also in the eddies on the edges of 
the layer. These eddies lie outside the high velocity 
gradient region of the layer as can be seen by compar- 
ing figure 18 with figure 22. Therefore, the resulting 
flame spreads transversely faster into the unreacted 
species than did the mixing layer defined by the high 
velocity gradient zone. Thus, the region of the mix- 
ing layer defined by the velocity gradient is not as 
transversely wide as the flame zone defined by tem- 
perature gradient in the mixing layer, in agreement 
with reference 16. 

Figures 23 through 29 show plots at seven stream- 
wise stations (x = 0.51, 0.58, 1.0, 2.0, 3.0, 4.0, 
and 5.0 cm) of the major chemical species (H 2 , O 2 , 
and H 2 O) and minor chemical species (OH, H, O, 
HO 2 , and H 2 O 2 ). Contour plots giving the two- 
dimensional distribution of the species are given in 
figures 30 through 35. Initially, at x = 0.51 cm 
(fig. 23(b)), fuel and air have just begun to mix, and 
no significant amount of water has yet formed. A 
very narrow band of hydrogen peroxide (H 2 O 2 ) is 
present just above the splitter plate center, and a 
very small amount of hydroperoxyl (HO 2 ) lies just 
below that spike. At x = 0.58 cm (fig. 24), the hydro- 
gen and oxygen profiles begin to broaden, but no wa- 
ter has yet appeared in the layer. The hydrogen per- 
oxide spike is still the most predominant, and, while 
the profile has not broadened, the peak has increased. 
(Note the ordinate in fig. 24(b) has been rescaled.) A 
small amount of hydroperoxyl still lies below the hy- 
drogen peroxide peak, and small amounts of atomic 
hydrogen (H) and atomic oxygen (O) have appeared 
there. At x — 1.0 cm, as described in figure 25, the 
hydrogen and oxygen profiles have developed, and a 
small amount of water (8 percent by mass) has been 
produced in a narrow profile below the splitter plate 
centerline. The H 2 and O 2 profiles are appropriately 
depressed in the region of the water peak. Notice- 
ably increased profiles of H, O, and OH also appear 
at this station just below the splitter plate centerline 
(y = 2.5 cm). The O and OH profiles lie slightly be- 
low the water peak, and the H profile lies just above 
the water peak; this is consistent with the spatial dis- 
tribution of reactant species. Small amounts of HO 2 
and H 2 O 2 still remain at and just above the plate 
centerline. 

Figure 26 diagrams the species profiles at x = 

2.0 cm. The H 2 , O 2 , and H 2 O profiles have broad- 
ened significantly more at this station, and the wa- 
ter peak has risen to approximately 23 percent by 
mass. The minor species profiles have also broad- 
ened significantly, with atomic oxygen peaking at 

3.0 percent, hydroxyl peaking at 2.0 percent, and 
atomic hydrogen peaking at 0.8 percent, all by mass. 
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Small amounts of hydroperoxyl and hydrogen per- 
oxide are still present just above the splitter plate 
centerline. With further movement downstream at 
x = 3.0 and 4.0 cm (figs. 27 and 28), the major and 
minor species profiles continue to develop, increasing 
both in width and in peak values. There are distinct 
distortions in the H 2 profiles in both figures because 
of eddies being located on the upper edge of the mix- 
ing layer. There is also a general migration of each 
profile to lower values of y with increasing streamwise 
coordinate. The increase in product species along 
the lower edge of the mixing layer is a direct re- 
sult of preferential burning in this region of the layer. 
The mixing layer is most nearly stoichiometric there, 
and the temperature reaches values that favor rapid 
ignition and combustion. At the last streamwise sta- 
tion given in figure 29, x = 5.0 cm, the major and 
minor species profiles broaden considerably further 
and shift to even smaller values of y . The noticeable 
increase in the rate of spread of the species profiles is 
associated with the second instability that is present 
in the mixing layer in this region and is consistent 
with transitioning to a turbulent state. 

Two-dimensional contour plots of the species are 
given in figures 30 through 35. The resulting struc- 
ture as the mixing layer develops, described previ- 
ously in figures 15 through 29, can be clearly seen 
in these figures. The first and second regions of in- 
stability are apparent for each species that is shown. 
The more rapid transverse spread of each species in 
the latter third of the layer can also be seen. A quies- 
cent region between the two instabilities also occurs 
for each species, as expected. Additionally, there is a 
general downward migration of each of the product 
species with increasing streamwise coordinate. The 
structure of the product species, typified by water, 
in the downstream region of the layer is also inter- 
esting. The vortical nature of the flow, seen earlier 
in figure 22, results in regions of unreacted hydrogen 
gas being captured by regions (or “folds”) of prod- 
uct water. Once captured, the regions of hydrogen 
have difficulty mixing with oxygen so that they can 
ultimately react. This phenomenon, often termed 
“unmixedness,” reduces the overall level of reaction 
that can be achieved, and contributes to a reduction 
in the efficiency of combustion. 

This completes the analysis of the spatially evolv- 
ing, supersonic, reacting mixing layer case using 
the finite-difference method. All the calculations 
described above were carried out on the Control 
Data Corporation VPS-32 computer (an expanded- 
memory CYBER 205) at the NASA Langley 
Research Center. The case required 5.1 CPU hours 
to reach the integration time of 0.1 m/s and used a 
core memory of 8 million 64-bit words. 


4.2 Simulations Using the Hybrid Spectral 
Algorithm 

The hybrid spectral algorithm has also been ap- 
plied to a spatially developing, laminar, supersonic, 
chemically reacting mixing layer. As noted earlier, 
no splitter plate dividing fuel and air is included in 
this case. Rather, initial profiles of flow variables are 
prescribed that approximate the flow some small dis- 
tance downstream of a splitter plate. Except for this 
modification, the configuration is identical to that 
considered in section 4.1. That configuration is de- 
scribed schematically in figure 36. Fuel is again intro- 
duced at 293 K, and air is introduced at 2000 K. Both 
fuel and air have an initial free-stream Mach number 
of 2.0, which ensures that no subsonic zone will occur 
in the mixing layer because of chemical heat addition 
or overall losses within the flow. The conditions re- 
sult in hydrogen and air velocities of 2604 m/s and 
1729 m/s, and Reynolds numbers of 4924 and 974, 
respectively. The previous study discussed in sec- 
tion 4.1 did contain a small subsonic zone in the im- 
mediate neighborhood of the splitter plate because 
of flow separation and a somewhat larger subsonic 
zone produced by heat addition in the later wake of 
the layer. It is advantageous to consider flows with 
the spectral method that are either fully supersonic 
or fully subsonic, as crossing a sonic line with the 
method requires special treatment. 

The overall domain considered in figure 36 is 
again 5 cm long, which allows sufficient length for 
initial development of the mixing layer. The do- 
main is mapped in the transverse direction to ±00 
with equation (111) of section 3.5.2. This ensures 
that the transverse boundaries lie well into the free 
stream so that the boundary conditions discussed in 
section 3.5.3 are properly posed. Initial and inflow 
boundary conditions are also chosen to be consistent 
with section 3.5.3, but in this case they are distrib- 
uted according to a hyperbolic tangent function that 
closely approximates the profiles of velocity, static 
temperature, and species that exit some small dis- 
tance downstream of a splitter plate. These profiles 
are also diagramed schematically in figure 36. 

By using the configuration and conditions de- 
scribed above, the resulting reacting mixing layer 
flow field was computed. Before detailed calcu- 
lations were begun, the hybrid spectral code was 
first checked against the earlier finite-difference code 
for this case. The calculations were carried out 
on a somewhat more coarse 51- by 51-point grid, 
and a simple one-step hydrogen-air chemistry model 
( 2 H 2 + O 2 2 H 2 O) was used in both codes to re- 
duce computation requirements. The detailed 18- 
equation chemistry scheme was common to both 
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programs, and therefore it did not require checkout 
in the spectral code. Results of the comparison at 
0.02 ms are given in figures 37 through 40. Agree- 
ment between the two programs is excellent in each 
plot, for both the fluid variables and species mass 
fractions. 

After the compatibility of the spectral and finite- 
difference codes was verified, the spectral code was 
then used to carry out more detailed calculations of 
the mixing layer flow of figure 36. Those calculations 
were performed on a grid of 201 points in the stream- 
wise (finite-difference) direction and 51 points in the 
transverse (spectral) direction. The grid was uniform 
in x and highly compressed in y in the region of the 
mixing layer. In fact, the spectral grid was chosen to 
be identical with the grid used for the finite- difference 
calculation in section 4.1 except for the streamwise 
compression employed about the splitter plate that 
was not included in the spectral calculation. The re- 
sulting flow field at 0.02 ms is described in figures 41 
through 53. 

Figure 41 shows a plot of streamwise velocity pro- 
files at four streamwise stations located at x = 0, 
1.0, 3.0, and 5.0 cm. Without the splitter plate to 
initiate disturbances and destablize the flow, there is 
only a small change in the profiles from the initial 
to the final streamwise station. Careful examination 
of figure 41 reveals the appearance of two-point os- 
cillations of small amplitude superimposed upon the 
velocity profiles. These are Gibbs oscillations that 
occur when a spectral method is used to resolve the 
large gradients that occur in this study. Gibbs os- 
cillations also occur when other numerical methods 
are employed, but the spectral method does not have 
sufficient numerical dissipation to damp the oscilla- 
tions. In this case, the numerical oscillations grow 
quite slowly with time, and a standard Laplacian fil- 
ter applied as a postprocessor following completion 
of the calculation is sufficient to remove them. Each 
dependent variable is therefore filtered by applying 
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where £ represents p, u, v, E , and /,*. When the 
coefficient leading the Laplacian is sufficiently small 
(Ay 2 /4 = 6 x 1CT 10 in this case) the filter dissipates 
only the two-point oscillations and leaves the values 
about their mean unaffected. Results from figure 41 
following filtering are given in figure 42. There are 
no structural changes in figure 42 relative to figure 41 
and the two-point oscillations apparent in figure 41 
have been removed in figure 42. 

Figure 43 shows a plot of streamwise velocity ver- 
sus streamwise coordinate at the three transverse sta- 


tions chosen in the finite-difference study. Note that 
the upstream and downstream instabilities that were 
present before do not occur in this case. The lack of 
the upstream instability is expected since the sepa- 
ration just downstream of the splitter plate initiated 
the instability in the finite-difference study. The lack 
of the downstream instability in this study seems to 
be due to the loss of a triggering mechanism for tran- 
sition by the upstream disturbance. To resolve this 
issue, the identical problem considered here was com- 
puted with the finite-difference code of section 4.1. 
Results from that analysis, again at 0.02 ms, are 
shown in figure 44. The results demonstrate the 
same behavior as those given in figure 43. There- 
fore, it appears that the downstream disturbance and 
transition of the mixing layer is dependent on an ini- 
tial triggering by the upstream instability. Signif- 
icantly longer calculations in time using the finite- 
difference program, in the absence of a splitter plate, 
never yielded transition within the 5-cm length of the 
physical domain of this problem. 

Figure 45 gives a plot of temperature versus 
streamwise coordinate at the three transverse sta- 
tions chosen in figure 43. Again, there is no evidence 
of the upstream and downstream instabilities present 
in the study with the splitter plate. Figure 46 de- 
scribes temperature profiles in the mixing layer at 
the four streamwise stations chosen in figure 42. As 
the flow evolves in x, there is a noticeable increase in 
temperature just below the center of the layer be- 
cause of the exothermic chemical reactions taking 
place there. The results of chemical reaction can be 
seen more directly in figures 47 through 53, which 
show the spatial evolution of the major (H 2 , O 2 , and 
H 2 O) and minor (H, O, OH, HO 2 , and H 2 O 2 ) chem- 
ical species at seven streamwise stations located at 
0, 0.4, 1.0, 2.0, 3.0, 4.0, and 5.0 cm. 

Figure 47 diagrams the initial distribution of re- 
actant species at x = 0 cm. No product species have 
formed at this station. Figure 48(a) shows the ma- 
jor species distribution a short distance downstream 
from the initial station at x = 0.4 cm. Initial re- 
action has begun at this location and a small nar- 
row profile of water (about 5 percent by mass) has 
formed. There are local depressions in the hydrogen 
and oxygen mass fraction profiles in the region of 
water production. Comparison of figure 48(a) with 
46 also shows that this region correctly corresponds 
to that of the peak temperature in the flow. As in 
the earlier finite-difference calculations of section 4.1, 
the zone of water production represents the region 
nearest stoichiometric conditions and at the required 
elevated temperature for chemical reaction, so it is 
reasonable that water initially forms here. The mi- 
nor species distributions at x = 0.4 cm are given 
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in figure 48(b). Atomic hydrogen lies in the imme- 
diate neighborhood of the water profile and peaks 
at 0.2 percent by mass, whereas the hydroxyl and 
atomic oxygen profiles extend well below the wa- 
ter and peak at 0.8 and 0.96 percent, respectively. 
Small amounts of hydroxyl (0.002 percent) and hy- 
drogen peroxide (< 10 -4 percent) also exist and lie 
just above the water peak. These distributions are 
again in agreement with the finite- difference results, 
with OH and O lying at or below the water where 
stoichiometry favors their higher population and H 
lying at or above the water for a like reason. 

Figure 49 shows major and minor species pro- 
files at 1.0 cm downstream from the initial sta- 
tion. All product species attain significantly higher 
peak values at this location, and the profiles have 
broadened considerably. Water peaks at 22 percent 
by mass, and hydroxyl, atomic hydrogen, and atomic 
oxygen peak at 2.4, 1.1, and 2.8 percent, respec- 
tively. The OH and O profiles shift to lower val- 
ues of y, whereas the H profile moves to a larger 
value. The hydroperoxyl profile shifts to a somewhat 
higher value of y but retains nearly the same peak 
value, whereas the hydrogen peroxide profile remains 
at nearly the same location and reaches a slightly 
higher peak (0.0007 percent). Species profiles con- 
tinue to broaden at x = 2.0 cm as shown in figure 50. 
There are slight increases in the peak values of wa- 
ter (23 percent) and atomic hydrogen (1.2 percent), 
and a slight decrease in atomic oxygen (2.6 percent). 
The remaining species remain essentially unchanged. 
There are very slight shifts in the transverse coordi- 
nate of the peak species values, but these shifts do 
not appear to be significant. 

There is no significant change in species profiles 
beyond the x = 2.0 cm station, and the chemistry 
appears to have reached a local equilibrium with the 
flow. Comparison of figure 50 at x = 2.0 cm with 
figures 51, 52, and 53 at x = 3.0, 4.0, and 5.0 cm, 
respectively, confirms this observation. This down- 
stream evolution of the chemistry differs considerably 
from the finite-difference study of section 4.1. This 
difference in evolution again appears to be linked 
to the absence of flow instabilities spawned by the 
presence of the splitter plate included in the finite- 
difference analysis. The effects of the splitter plate 
on flow instabilities were discussed in section 3 and in 
section 4.1. Without the presence of the first insta- 
bility, early mixing is reduced and transition does not 
occur in the layer within the limits of the domain that 
is considered. In the absence of downstream transi- 
tion, downstream mixing is significantly retarded. 

To study the effects of the inflow perturbations 
imposed on the flow by the splitter plate, the finite- 
difference calculation of section 4.1 was reconsidered. 


It was found from analyzing computed results as a 
function of time at the first station downstream of 
the splitter plate trailing edge that perturbations 
imposed on the flow variables could be correlated 
quite well by using a single perturbation function. 
That function was given by 

= e~( ay ^ 2 A sin (wt) (115) 

where 

a = 1000 
A = 0.064 

oj = 12 271061 rad/s 

The exponential term damps the effects of the per- 
turbation with transverse movement away from the 
plate. The trigonometric function describes the pe- 
riodic nature of the disturbance, and A is the am- 
plitude of the disturbance, normalized by the free- 
stream velocity. 

The flow perturbation described by equa- 
tion (115) is now applied to the present analysis by 
using the spectral program. Recalling the discussion 
in section 1, it was noted in reference 21 that for 
reacting flows, the eigenvalues that determine flow 
stability are only weakly dependent on Mach num- 
ber. Therefore, it appears reasonable to apply the 
perturbation data from the finite-difference analysis 
to this problem. Each primitive fluid variable was 
then perturbed at the inflow boundary as described 
by equation (115). The initial streamwise velocity, 
for example, is then given by 

u = £700(1 + 0') (116) 

A similar procedure is also used to describe the 
inflow density and pressure. All other required fluid 
inflow data can then be computed as usual. The 
transverse velocity component remains fixed at zero 
as in the previous analysis, roughly representing the 
flow just downstream of the recompression shock in 
the finite-difference study. The spectral technique 
cannot capture strong shocks in supersonic flow that 
would occur with the imposition of the transverse 
velocity, and, therefore, that problem could not be 
dealt with here. 

Results from the spectral study, again at 0.02 ms, 
using the inflow perturbation described above are 
given in figures 54 through 70. Figure 54 shows a 
plot of streamwise velocity versus streamwise coordi- 
nate at the same three transverse stations pictured 
in figure 43 for the unperturbed study. Note that the 
instability introduced at the inflow boundary now 
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persists through the solution domain. The distur- 
bance does not appear to be amplified significantly, 
however. There are three regions of instability in the 
streamwise direction, each separated by a zone where 
the oscillations are damped. These damped zones oc- 
cur at approximately 0.7 and 3.0 cm. These results 
are similar in certain respects to those seen in the 
finite-difference analysis at 0.02 ms as given in fig- 
ure 16(d). There are three regions of instability in fig- 
ure 16(d), but the second damped zone occurs further 
upstream. Also, the downstream region of instabil- 
ity is larger, of greater amplitude, and of increasing 
wavelength. From this comparison, it appears likely 
that the mechanism initiating the instabilities in the 
finite-difference study is more complex than that as- 
sumed in the present study. The assumed form of 
the perturbation does, however, allow study of the 
effects of such an instability on the development of 
the mixing layer and the resulting chemical reactions. 
These effects are assessed by comparing results from 
perturbed and unperturbed spectral studies. 

Figure 55 gives a plot of streamwise velocity ver- 
sus transverse coordinate at four streamwise stations. 
The unperturbed results at the same four stations are 
given in figure 42. Note that now there is a marked 
shift in the overall magnitude of each velocity profile 
because of the temporal perturbation, but the profile 
development is not significantly affected. 

Figure 56 shows a plot of streamwise velocity con- 
tours in the mixing layer. The structure described 
by figure 54 is shown to persist throughout the layer, 
and there is no marked growth of the mixing layer 
thickness as defined by the velocity gradient. Sim- 
ilar results are yielded by the plot of vorticity con- 
tours given in figure 57. There again is no signifi- 
cant growth of the layer, but the perturbation does 
produce vortical structure, albeit of lesser scale than 
that observed in the calculation including the splitter 
plate. 

Figure 58 describes the temperature field with in- 
creasing streamwise coordinate at the three trans- 
verse locations used in the unperturbed solution of 
figure 45. The streamwise structure is again appar- 
ent in this plot. Comparison of figure 58 with fig- 
ure 19, which gives the streamwise temperature de- 
velopment in the flow with the splitter plate, shows 
similar behavior of the profiles near, just above, and 
below the plate. The profiles are quite different well 
downstream of the plate, however, because of the ab- 
sence of transition without the plate. Figure 59 de- 
scribes transverse temperature profiles with stream- 
wise distance in the perturbed layer. The results are 
quite similar to those of the unperturbed flow given 
in figure 46, although slightly higher peak tempera- 
tures are achieved in the perturbed solution. 


Profiles and contours of the chemical species 
present in the perturbed reacting mixing layer are 
given in figures 60 through 70. Figure 60 shows pro- 
files of the major species (H 2 , O 2 , and H 2 O) at the 
initial station of the calculation. No water has been 
formed at this station. Since the chemical species 
are unperturbed at the inflow boundary, figure 60 
is identical to the results for the nonperturbed case 
given in figure 47. Figure 61 describes the species 
distribution at x = 0.4 cm downstream of the in- 
flow boundary. The water profile is somewhat more 
broad but has a slightly lower peak than the unper- 
turbed solution in figure 48(a). The minor species 
also peak about 0.1 percent lower than the unper- 
turbed solution given in figure 48(b). The small dif- 
ferences in the two solutions appear to be due to a 
somewhat higher rate of mixing due to the perturba- 
tion, which would increase the transverse spread of 
the profile and reduce its peak. The trends in pro- 
file spread established at x — 0.4 cm continue at the 
x = 1.0-cm station given in figure 62. By compar- 
ing figure 62 with 50, it can be seen that the profile 
peaks are identical, but the perturbed profiles are 
slightly broader. With each succeeding downstream 
station beyond x — 1.0 cm, the species profiles of the 
perturbed problem continue to become transversely 
more broad than the unperturbed study because of 
the improved mixing afforded by the perturbation. 
The mixing process remains laminar, however, since 
the induced instability is never sufficient to trigger 
transition in the latter portion of the mixing layer. 
This behavior can be seen even more clearly in the 
product species contour plots of the perturbed mix- 
ing layer, given in figures 67 through 70. The per- 
turbation on the fluid variables induces an instabil- 
ity in these species that is initially quite strong. The 
instability decays with downstream movement, how- 
ever, and it has essentially dissipated by the time the 
outflow boundary is reached. Therefore, the species 
distributions indicate, as did the fluid variable re- 
sults, that the mechanism that triggers transition in 
the mixing layer flow is more complex than the func- 
tion that was assumed. These issues of reacting flow 
stability are further addressed in the conclusions dis- 
cussed in section 5. 

The calculations described above were carried 
out on the VPS-32 computer at the NASA Langley 
Research Center. The case required 7.0 CPU hours 
to reach the integration time 0.02 ms and used a 
core memory of 7.2 million 64-bit words on a 201 
by 51 grid. 

5. Conclusions 

Research has been undertaken in this study to 
achieve an improved understanding of important 
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physical phenomena present when a supersonic flow 
undergoes chemical reaction. To explore the behavior 
of such flows, detailed physical models of convective 
and diffusive mixing and finite-rate chemical reaction 
in supersonic flow were developed. Two numerical 
algorithms were then constructed to solve the equa- 
tions governing supersonic chemically reacting flow 
that resulted from these models. The first algorithm 
was developed around an established finite-difference 
technique modified to consider multicomponent re- 
acting flow. The second algorithm employed a hybrid 
pseudospectral technique in one spatial direction for 
improved resolution of the reacting flow field. The 
previous scheme was retained in the other spatial di- 
rection. Computer programs were written using both 
algorithms, and each program was used to study a 
spatially developing and reacting supersonic mixing 
layer. The results obtained from these studies were 
then analyzed, and conclusions were drawn concern- 
ing the structure of the reacting mixing layer. Those 
conclusions, which were discussed in section 4, are 
now summarized. 

Supersonic reacting flows exhibited many of the 
same features observed for subsonic reacting and 
nonreacting flows. In particular, the vortical struc- 
ture of the flow, noted in much of the subsonic nonre- 
acting flow literature, was shown for the first time to 
be quite predominant in supersonic reacting flow as 
well. In agreement with the earlier reacting subsonic 
literature, the vortical structure had a marked effect 
on chemical reaction in supersonic flow. Significant 
burning took place in the eddies on the edges of the 
mixing layer, broadening the reaction zone relative to 
the layer thickness defined by the velocity gradient. 
In addition, the vortical flow resulted in the roll up of 
unburned reactants inside a layer of partially or fully 
burned products. This phenomenon, often termed 
“unmixedness” in subsonic flows, prohibited the re- 
action of captured reactants and reduced the overall 
efficiency of the combustion process. Unmixedness 
was thus shown for the first time to be a potential 
problem in reducing the efficiency of supersonic com- 
bustion as well as subsonic combustion, and tech- 
niques will likely be needed to reduce its effects. 

Calculations with the present model also showed 
that at supersonic speeds the reacting mixing layer 
remained laminar for the region studied if no external 
disturbance to trigger transition to turbulence was 
introduced. When a splitter plate was used to ini- 
tially separate fuel and air, however, it provided the 
required disturbance. The unstable recirculating flow 
that formed at the base of the splitter plate, follow- 
ing the Prandtl-Meyer expansion off the plate and the 
unstable recompression shock a short distance down- 


stream, provided that disturbance. The resulting os- 
cillatory flow then propagated downstream triggering 
transition-like phenomena in the latter fourth of the 
domain being studied. Mixing of fuel and air then 
improved dramatically in this region, markedly in- 
creasing chemical reaction as evidenced by the spread 
of product profiles. To study the effect of heat re- 
lease in this region, calculations were also carried 
out without chemical reaction. Results showed that 
the unstable region near the splitter plate was un- 
affected when reaction was removed. There was no 
reaction in the early part of this region, and reac- 
tion was mainly endothermic further downstream, 
so little effect was expected. Well downstream in 
the transition-like region, the reaction was highly 
exothermic, however, and removing chemical reac- 
tion (and therefore, chemical heat release) caused the 
amplitude of the disturbance there to increase signif- 
icantly. This result was in agreement with earlier 
experimental and numerical literature for subsonic 
flow, where it was found that heat release retarded 
mixing. This effect was thus shown for the first time 
to occur in supersonic reacting flow as well. 

This study also represented the first application of 
spectral methods to study supersonic reacting flows. 
The hybrid spectral method employed in this study 
was used to predict the spatial development of a su- 
personic, chemically reacting mixing layer. The first 
case studied considered the development of a mix- 
ing layer downstream of a splitter plate separating 
fuel and air. No plate was included in the calcu- 
lation; rather the effects of the plate were modeled 
by using an appropriate initial profile. As in one 
of the finite-difference studies, the layer without the 
plate never developed a sufficient upstream distur- 
bance to trigger transition in the downstream region 
of the problem that was studied. To initiate transi- 
tion, data were taken from the upstream disturbance 
that caused transition in the finite-difference study 
and were correlated to form an initial perturbation 
function on the inflow field of the spectral study. The 
perturbation alone was not sufficient to trigger tran- 
sition in the spectral study, although species mix- 
ing and chemical reaction were enhanced well down- 
stream. It was therefore concluded that transition 
was initiated in the finite-difference study by a mech- 
anism more complicated than that represented by 
the simple perturbation function used in the spec- 
tral study. 


NASA Langley Research Center 
Hampton, VA 23665-5225 
September 13, 1988 
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Table I. Finite-Rate Chemistry Model and Rate Coefficients for Each Reaction 


Reaction 

number 

Reaction 

Reaction rate variables 

Ai 

Ni 

Activation 

energy, 

cal/g-mole 

1 

H 2 + 0 2 = OH + OH 

.170 x 10 14 

0 

48150 

2 

H + 0 2 = OH + 0 

.142 x 10 15 

0 

16 400 

3 

OH + H 2 = H 2 0 + H 

.316 x 10 8 

1.80 

3 030 

4 

O + H 2 = OH + H 

.207 x 10 15 

0 

13 750 

5 

OH + OH = H 2 0 + O 

.550 x 10 14 

0 

7000 

6 

H + OH = H 2 0 + M 

.221 x 10 23 

-2.00 

0 

7 

H + H = H 2 + M 

.653 x 10 18 

-1.00 

0 

8 

H + 0 2 = H0 2 + M 

.320 x 10 19 

-1.00 

0 

9 

ho 2 + oh = h 2 o + o 2 

.500 x 10 14 

0 

1000 

10 

H0 2 + H = H 2 + 0 2 

.253 x 10 14 

0 

700 

11 

H0 2 + H = OH + OH 

.199 x 10 15 

0 

1800 

12 

H0 2 + O = OH + 0 2 

.500 x 10 14 

0 

1000 

13 

ho 2 + ho 2 = H 2 0 2 + 0 2 

.199 x 10 13 

0 

0 

14 
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Figure 1. Cross section of transitional “vortex” eddy in gas diffusion flame (ref. 11). 
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Figure 2. Cross section of large turbulent eddy in gas diffusion flame (ref. 11). 





Figure 5. History of oxygen mass fraction. 



Figure 6. Histories of hydroxyl and water mass fractions. 
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Figure 7. Axial velocity distributions. 



0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 

x, m 

Figure 8. Axial temperature distributions. 
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Figure 9. Axial pressure distributions. 
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Figure 11. Axial oxygen mass fraction distributions. 
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Figure 12. Axial hydroxyl and water mass fraction distributions. 
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Figure 13. Steady-state residual reduction rates of the methods. 
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Figure 15. Mixing layer velocity vector field (every fourth vector plotted (0.1 ms)). 
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(a) 0.1 ms. 

Figure 16. Streamwise velocity versus x at y locations. 
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(b) 0.09 ms. 
Figure 16. Continued. 
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(d) 0.02 ms. 
Figure 16 . Continued. 










Figure 17. Streamwise velocity versus y at x locations. 



Figure 18. Streamwise velocity contours in mixing layer. 
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Figure 21. Temperature contours in mixing layer. 



Figure 22. Vorticity contours in mixing layer. 
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(a) Major species. 

Figure 23. Mass fraction versus y at z — 0.51 cm. 
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(b) Minor species. 
Figure 23. Concluded. 
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(a) Major species. 

Figure 25. Mass fraction versus y at x — 1.0 cm. 
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(a) Major species. 

Figure 26. Mass fraction versus y at x = 2.0 cm. 
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(b) Minor species. 
Figure 26 . Concluded. 
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(b) Minor species. 
Figure 27. Concluded. 
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(a) Major species. 

Figure 28. Mass fraction versus y at x = 4.0 cm. 
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(b) Minor species. 
Figure 28. Concluded. 
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(a) Major species. 

Figure 29. Mass fraction versus y at x — 5.0 cm. 





Figure 30. Hydrogen mass fraction contours in mixing layer. 



Figure 31. Oxygen mass fraction contours in mixing layer. 



Figure 32. Atomic hydrogen (H) mass fraction contours in mixing layer. 



Figure 33. Atomic oxygen (O) mass fraction contours in mixing layer. 



Figure 35. Water mass fraction contours in mixing layer. 



Figure 36. Reacting mixing layer schematic for spectral calculations. 
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(a) Major species. 

Figure 48. Mass fraction versus y at x = 0.4 cm. 
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(b) Minor species. 
Figure 48. Concluded. 
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(a) Major species. 

Figure 49. Mass fraction versus y at x = 1.0 cm. 
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(a) Major species. 

Figure 50. Mass fraction versus y at x = 2.0 cm. 
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(a) Major species. 

Figure 51. Mass fraction versus y at x = 3.0 cm. 











(a) Major species. 

Figure 53. Mass fraction versus y at x = 5.0 cm. 
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(b) Minor species. 
Figure 53. Concluded. 
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Figure 54. Streamwise velocity versus x at y locations following inflow perturbation. 
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Figure 55. Streamwise velocity versus y at x locations. 






Figure 56. Streamwise velocity contours in mixing layer. 



Figure 57. Vorticity contours in mixing layer. 
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(a) Major species. 

Figure 62. Mass fraction versus y at x — 1.0 cm 
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(a) Major species. 

Figure 63. Mass fraction versus y at x = 2.0 cm. 
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(a) Major species. 

Figure 65. Mass fraction versus y at x = 4.0 cm. 






fraction 
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Figure 66. Mass fraction versus y at x = 5.0 cm. 





Figure 69. Hydroxyl mass fraction contours in mixing layer. 



Figure 70. Water mass fraction contours in mixing layer. 
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